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ABSTRACT 


Wider acceptance of CFD methods in the design and the analysis require development of 
CFD techniques that should satisfy two main criteria; provide results within a reasonable time 
frame and with apriori known accuracy. Major efforts have been spent in recent years to establish 
the accuracy of numerical solvers. However more research is required to develop reliable and 
robust numerical solvers suitable for modeling of complex turbomachinery flows. 

An unsteady viscous flow solver based on the Runge-Kutta scheme has been developed. 
Pseudo-time step technique has been incorporated to provide efficient simulation of unsteady 
flows. Utilization of the pseudo-time approach reduces the computational time by a factor varying 
from 5 to 25 times in comparison with the original solver. The results of the stability analysis of 
the dual time step scheme are used to establish an optimum pseudo-tune step based on the local 
CFL, VonNeuman numbers and the ratio of pseudo-to-physical time steps. Code has been modified 
to incorporate multi-block capabilities. This modification is essential for the modeling of complex 
multidomain configurations such as the rotor-stator interaction, film cooling, etc. Multoblock 
feature simplifies the complexity of the grid generation process. It also improves the quality of the 
grid, thus contributing to the enhanced accuracy of the numerical modeling. 

The code has been validated against the analytical and experimental data. The influence of 
the numerical aspects (artificial dissipation, grid density etc.) on the accuracy of the prediction of 
the wake decay, transition, flow over a cylinder has been analyzed. Based on the results of the test 
cases, modifications to the k-e model to improve the accuracy in the regions with dominant normal 


stresses are incorporated. 
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Numerical simulation of the unsteady flow in a compressor cascade has been performed to 
assess the ability of the code to simulate unsteady flow in a turbomachinery cascade. Three low Re 
k-e turbulence models have been assessed for their ability to predict the unsteady transitional 
flows. Good agreement with the measured data and with an earlier Euler/boundary layer prediction 
has been achieved. The numerical solver was able to predict major features, associated with the 
wake induced transition on a compressor blade (wake induced transitional strip, wake induced 
turbulent strip, etc.). Analysis and interpretation of the results from the unsteady flow simulation 
have been carried out to understand additional flow physics associated with the transitional flow. 

A coupled experimental and computational study of the effects of the nozzle wake-rotor 
interaction in a turbine is carried out to understand the cause and effects of the unsteady flow in 
turbine rotors. The result of the numerical prediction correlates well with the Laser Doppler 
Velocimeter data and dynamic pressure measurements. An assessment of the viscous and the 
inviscid contribution to the nozzle wake decay and the unsteady loss distribution in the rotor 
passage reveals the dominant effect of the viscous decay upstream of the leading edge. Inside the 
passage, the inviscid effects have a significant influence. The predicted flow at the off-design 
condition has been interpreted to understand the nature of the unsteady flow field at a high negative 
incidence angle. 

Variation of the flow Reynolds number between the take off and cruise conditions 
significantly affects the boundary layer development on a low-pressure turbine blading. A 
decreased Reynolds number leads to flow separation on the suction surface of the blading, thus 
increasing losses. A numerical simulation has been carried out to assess the ability of the Navier- 
Stokes solver to predict transitional flows in a wide range of Reynolds numbers and inlet 
turbulence intensities. A number of turbulence (including the Algebraic Reynolds Stress Model) 


V 


and transition models have been employed to analyze the reliability and accuracy of the numerical 
simulation. A comparison between the prediction and the experimental data reveals a good 
correlation. However, the analysis shows that the artificial dissipation in the numerical solver may 
have a profound effect on the transition in a separated flow. 

The viscous flow solver has been employed for the numerical investigation of the 
aerothermal field due to the leading edge film cooling at a compound angle. Good agreement with 
the measured data has been achieved. Results of the numerical investigation have been used to 
analyze the vortex- structure associated with the coolant jet-freestream interaction to understand the 
effect of different vortices on the cooling effectiveness and aerothermal losses. Two counter- 
rotating vortices generated by the interaction between the mainflow and the coolant jet have been 
found to have a major influence in decreasing the cooling efficiency through strong entrainment of 
the hot fluid. Results of the numerical simulation indicate that the turbulence length scale has a 
significant effect on the accuracy of the numerical prediction of the film cooling. Not only the inlet 
turbulence intensity but also the turbulence length scale should be accurately prescribed to achieve 
reliable numerical prediction of the heat and the mass transfer due to the film cooling. 

Numerical analysis of the tip leakage flow in a turbine is utilized to investigate physics of 
the seco ndar y flow in a rotor including interaction between secondary vortices, tip leakage vortex, 
and rotor wake. Analysis of the leakage flow development shows that the relative motion of the 
blade and the casing wall reduces the propagation of the leakage flow into the mainflow. The tip 
leakage vortex is confined to the suction surface comer of the casing. Most of the leakage losses is 
due to the mixing of the tip leakage vortex downstream of the trailing edge. 
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Chapter 1 


INTRODUCTION 


1.1 Significance 

Significant research efforts and design advances have led to thermal efficiency of 
up to 60% in modem gas turbines. The modem compressor stage has efficiency of about 
90% and the modem turbine stage has efficiency of up to 95%. Further improvements in 
efficiency become more and more difficult and require a much deeper understanding of the 
flow field inside turbomachines. 

Three methods have been widely used in the analysis and design of modem 
turbomachinery; experimental research, analysis in conjunction with empirical data base, 
and numerical methods. The development of the accurate numerical methods and 
computer hardware have made the numerical simulation more efficient, reliable and 
affordable. However, questions of affordability and reliability of the numerical simulation 
are the key factors in the future development of the numerical analysis for the 
turbomachinery design and analysis. 

Wider acceptance of CFD methods in the design and the analysis requires CFD 
techniques that should satisfy two main criteria; provide results within reasonable time 
frame and with apriori knowledge of accuracy. Ideally, engineer need to have a technique 
(i.e. implementation of correct physical models etc.) and solution with a fixed level of 
precision. Despite complexity of the process, it is possible to establish reliability of the 
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numerical technique. Reliability and suitability of the mathematical model, including such 
aspects as the turbulence model, as well as mutual influence of the model and the 
numerical technique can not be verified for a general case. Major efforts have been spent 
in recent years to establish the accuracy of numerical models utilizing the benchmark data. 
These attempts had only limited success. Further development of the CFD requires more 
systematic approach to the assessment of the numerical accuracy. Melnik et al., 1995, 
suggested three major steps in assessing the capability of the CFD code; code verification, 
validation and certification. Code verification against known analytical solutions, reveals 
numerical accuracy of the numerical technique used. The second step is the code 
validation. This step is critical in assessing the ability of the code to provide an accurate 
solution for the benchmark test cases in a relatively narrow range of flow features. Last 
step is the code certification, which includes predictive capabilities of the code for 
complex and realistic cases. This step includes a systematic simulation of flow cases and a 
comparison with the existing experimental data. 

Flows in turbomachinery blade rows are very complex. They are strongly three- 
dimensional, viscous, with several types of secondary flows and vortices (horseshoe 
vortex, passage vortex, leakage flow, etc.). Their interaction with the blade, boundary 
layers and wakes results in mixing losses. Transitional flows and the high turbulence 
intensity result in additional complexities. Because of the complex nature of the flow, 
analytical methods are scarce and not accurate. Results of experimental investigations are 
limited to a narrow range of flow parameters in modeled turbomachines. 
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One of the major features of the flow in the blade passage is the unsteadiness 
caused by the relative motion of the stator and rotor. While the unsteadiness plays an 
important role in the flow through turbomachinery blade rows, majority of flow 
simulations are carried out assuming a steady flow approximation. The inlet flow pattern is 
prescribed as uniform and steady. This approach leads to the neglection of such 
phenomena as a rotor/stator interaction, vortex shedding and other unsteady flow effects. 
The unsteadiness has a major influence on the surface pressure distribution and shear 
stresses at the wall. The unsteady dynamic and thermal loading can reduce the life of the 
blades. To ensure reliable operation, the natural frequency of the blade should be different 
from the frequency of the vibration caused by the flow unsteadiness. The main sources of 
the unsteadiness are; potential effect, wake-blade interaction, vortex-blade interaction 
random unsteadiness of the mean flow. 

Another problem which is closely connected to the unsteady nature of the flow in 
turbomachines is the heat transfer. Prediction of the heat transfer and film cooling effect is 
crucial to an understanding of the turbine flow field. Excessive blade temperature may lead 
to a thermal fatigue. Accurate analysis of this phenomenon is essential for good design. 

1.1.1 Unsteady Flow 

Unsteady interaction increases losses, blade vibration, and noise generation; and 
affects heat transfer in turbines. An understanding of the physics of the unsteady flow will 
enable an improvement in the overall aerodynamic and mechanical performance of the 
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turbine. Increases in the available computer resources and the development of more 
efficient computational algorithms in recent years have made the numerical simulation of 
unsteady flows more affordable. 

There are three main approaches to unsteady calculations in the blade passage. The 
utilization of the linearized Euler equation is the earliest and the least numerically intensive 
approach. The simulation based on the full Euler equation or the coupled Euler/boundary 
layer approximation, is used in the second approach. The last approach is the numerical 
simulation using Reynolds averaged Navier-Stokes equations with an adequate turbulence 
model. 

The linearized inviscid theory is based on the approximation that the perturbations 
of the mean flow are small and hence, these parameters are presented through Taylor 
series neglecting higher order terms. The more complex approach is the modeling of the 
flow using non-linear Euler equations. In this method, fall Euler equation is used and only 
the inviscid approximation is invoked. This approach was used in Giles (1988), He (1989), 
Domey and Verdon (1994), Fan and Lakshminarayana (1994). Giles (1988) analyzed the 
flow using the Lax-Vendroff scheme with non-reflecting boundary condition. He (1989) 
used a multi-step Runge-Kutta scheme to simulate two dimensional flow over an 
oscillating blade. A moving grid zone was implemented near the blade surface to treat 
blade oscillations. One of the approaches which combines the advantages of both the Euler 
and the boundary layer method is due to Fan and Lakshminarayana (1994). Fan and 
Lakshminarayana (1994) used a multi-step Runge-Kutta scheme with non-reflecting 
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boundary conditions. The results of the inviscid solution were used as an input for 
unsteady boundary layer calculations. Numerical results showed a good agreement with 
the experimental data. 

The last and the most complex approach is the Navier-Stokes simulation. Codes 
based on this approach were used by Rai (1987), Domey and Davis (1991), Ho and 
Lakshminarayana (1993), Fourmaux (1994), Amone et al. (1994), One of the earliest 
works in this field is the simulation done by Rai (1987). He used a thin layer 
approximation of the Navier-Stokes equations and a third order upwind difference 
scheme. An O-type overlaid grid was used. Ho and Lakshminarayana (1993) developed an 
unsteady Navier-Stokes solver based on a pressure correction method. Code was validated 
for the grid sensitivity and artificial dissipation. Numerical simulation of the rotor-stator 
interaction showed good correlation with experimental results. Fourmaux (1994) 
implemented four-step Runge-Kutta numerical scheme with combined H-O type grid. 
Amone et al. (1994) also applied the explicit four-step Runge-Kutta scheme where, for 
economy, the viscous terms were evaluated only on the first stage. 

A number of additional problems must be solved when the numerical modeling is 
based on unsteady Navier-Stokes equations. Unsteady numerical simulation results in a 
significant increase in required CPU time. Wave dissipation and dispersion characteristics 
of the steady state numerical scheme are not suitable for unsteady flow simulations. 
Problems associated with efficiency and accuracy of unsteady simulations can be 
overcome through the utilization of dual-step approach. The utilization of the inner cycle 
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to relax a time step limitation of the explicit schemes (e.g. Amone, 1994), or to remove 
the linearization error in the implicit scheme (Rai, 1987) was found to be essential for 
unsteady simulations. 

The choice of physically realistic turbulence model suitable for the unsteady 
calculation is a difficult one. Most of the turbulence models have been developed for the 
steady flow. Special attention needs to be paid to the ability of the model to resolve time 
scales associated with the flow unsteadiness. Many authors used simple models which had 
been validated only for steady flows, such as the 2 layer Baldwin-Lomax turbulence model 
(Amone et. al. (1994) and Domey et al. (1994)). Partially, this may be attributed to the 
additional stability and convergence problems caused by the incorporation of more 
complex turbulence models in unsteady computations. Turbulence models, as well as their 
improvements for unsteady flows, need to be more carefully investigated. Fan, 
Lakshminarayana, and Barnett (1993) modified two-equation k-e model for a application 
in the unsteady flow and showed good agreement with the experimental data for the flat 
plate and cascade unsteady flows. 

Even though many attempts have been made to develop and use unsteady Navier- 
Stokes solvers, none have been satisfactory validated against accuracy, especially in 
respect to the unsteady viscous layer near blade and wall surfaces. This is one of the major 
objectives of this thesis. 
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1.1.2 Transition to turbulence 

One of the challenging problems in turbomachinery is to understand the flow 
physics due to transitional flows associated with the laminar separation and the rotor- 
stator interaction in low pressure turbines. The rotor-stator interaction flow is inherently 
unsteady and transitional. Additional complexities arise due to these transitional boundary 
layers along the blade surfaces. Such complex unsteady and transitional boundary layer 
flow is known to affect the aerodynamic and thermal performance of a turbomachine 
(Simon and Ashpis, 1996). The transition from the laminar to turbulent flow on the blade 
surface is a common, yet complex, phenomenon in turbomachinery. The boundary layer 
development, losses, efficiency, and heat transfer are greatly affected by the transition. The 
ability to accurately predict the onset and length of the transition is very important in the 
design of efficient machines. 

There are three types of transition. The first is called the “natural” transition, 
where the laminar boundary layer develops the Tollmein-Schlichting wave, followed by an 
amplification of instabilities and finally the fully turbulent flow. Natural transition usually 
occurs with small ffeestream disturbances. The second type of transition is caused by large 
external flow disturbances. It is called the “bypass” mode because there is no Tollmein- 
Schlichting instability. The third type is called the “separated-flow” transition, which 
happens within the laminar boundary layer separation and may or may not involve the T-S 
wave. In turbine flows, the ffeestream turbulence level is usually high. Transition in these 
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flows is of the bypass mode most of the time. Natural transition is almost non-existent in 
practical flows. Separated flow transition is also common in low pressure turbines. 

The periodic passing of upstream wakes can also lead to transition patches on the 
downstream blade surface. This modification of the transition process is called the “wake- 
induced” transition. When the upstream wake impinges on the downstream blade surface, 
within a laminar boundary layer, transition occurs because of the sudden and large 
disturbance caused by the wake and the high turbulence level inside the wake. Periodic 
turbulence and transition patches may develop and transport downstream at the certain 
fraction of the period while laminar regions exist at the rest of the period. 

High performance and durability of turbines can be realized through an improved 
understanding of the physics associated with the transition and rotor-stator interaction 
phenomena. Considerable attention has been paid in recent years in developing computer 
codes to predict unsteady aerodynamics and heat transfer, but these efforts are hampered 
by a lack of understanding of the basic physics associated with these interactions and the 
lack of adequate physical modeling (transition/turbulence models), and validation of the 
codes. The ultimate solution of this problem has to come from a systematic, scientific, and 
building block approach. The measurement in an actual engine is not only complicated, but 
will rarely provide an insight into numerous sources of unsteadiness and mechanisms. 
Likewise, a computational code with the artificial dissipation and the numerical error may 
mask some of the important physics. The code has to be validated at several stages to 
ensure that the flow physics is captured accurately. The past computational effort was 
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main ly concerned with the large code development for the steady and the unsteady viscous 
flow in turbomachines using simple algebraic eddy viscosity models. Even the steady flow 
prediction with higher order turbulence models is not satisfactory due to inadequacy of 
physical models employed. The models do not adequately account for effects of rotation, 
curvature, heat transfer, compressibility, three-dimensional strain field, flow separation, 
and the unsteady flow. These codes are among the most sophisticated and comprehensive 
available for turbomachinery flows. There is a need to assess these numerical techniques 
and improve the computational efficiency. 

Mayle (1991) reviewed the transition phenomena in gas turbines. From a 
theoretical perspective, first transition is viewed as a sudden jump from the laminar to the 
turbulent flow. Laminar flow is separated from the turbulent flow by a single line or 
section. Through the modification to the boundary layer properties, two zones are patched 
together. This is the approach most numerical methods adopt (Launder and Spalding, 
1974; Schmidt and Patankar, 1991). 

Experimental results show that that the transition is not an abrupt process. 
Emmons (1951) is the first to propose that the transition is a three-dimensional and 
unsteady process, which has a region where laminar and turbulent flows co-exist. At a 
certain point in space, the flow could be laminar at sometime and turbulent at other times. 
This is the “intermittency” phenomenon. Most of the earlier theoretical investigations were 
concerned with the intermittency factor. Narasimha (1957) modified the Emmons theory 
through the hypothesis of concentrated breakdown and showed good agreement between 
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his model and the intermittency measured by Schubauer and Klebanoff (1955). Dhawan 
and Narasimha (1958) developed a model for the intermittancy based on a Guassian 
distribution. This model can only be used in combination with analysis/computational 
techniques to predict onset location and transition length. Mayle and Paxson (1991) later 
proposed a new theory that accounts for the extra terms due to the interaction between 
non-turbulent and turbulent flows in molecular and heat flux stresses. Models based on the 
experimental data are used to derive onset and transition. These models are formulated to 
account for the effects of turbulence intensity (Gostelow and Blunden, 1989; Abu- 
Ghannam and Shaw, 1980), pressure gradient (Gostelow and Walker, 1991), and other 
factors like curvature and surface roughness. 

Computations of transitional flows can be classified into four groups. The simplest 
one is a linear combination model. Transitional flow is divided into turbulent and laminar 
parts according to the intermittency parameter. Predictions using the model by Dhawan 
and Narasimha (1958) give an excellent agreement for two-dimensional flows without 
pressure gradient. But the extension to more complex flow situations has not been 
successful. The second method is incorporated in the framework of algebraic turbulence 
models. The total viscosity is assumed to be the sum of the molecular viscosity and the 
product of the intermittency and the eddy viscosity. The third group employs the one- or 
two-equation turbulence models, and will be discussed in the next paragraph. The fourth 
method uses the direct numerical simulation of three-dimensional unsteady flows. No 


models are needed in this case. 


Computation of the flow field including the transition by two-equation turbulence 
models is a popular approach (Jones and Launder, 1972; Schmidt and Patankar, 1991; Fan 
and Lakshminarayana, 1996). The low-Reynolds-number form of two-equation models are 
capable of capturing the transition inception location to a certain accuracy. Schmidt and 
Patankar (1991) examine the effects of inlet locations, inlet profile of the turbulence 
kinetic energy and the dissipation rate, and the freestream turbulence intensity on the 
transition using low-Reynolds-number two-equation models by Jones and Launder (1972), 
and Lam and Bremhost (1981). Through the modification of the turbulence production 
term and the introduction of two parameters to their model, onset and end of transition on 
a flat plate are accurately predicted. Fan and Lakshminarayana (1996) proposed a new 
model which modifies near wall functions and obtaine improved wake-induced transitions. 

There has been very limited computational effort to resolve the flow physics and 
provide an accurate prediction of the unsteady viscous layers in turbines. Fan and 
Lakshminarayana (1996) used an Euler-boundary layer approach and modified the 
turbulence models to predict unsteady transitional viscous layers in compressors and 
turbines for which detailed data is available (Schulz et al., 1990 and Halstead et al., 1995). 
No attempts have been made to assess the capability of the Navier-Stokes code to predict 
these unsteady transitional viscous layers due to the wake-blade interaction. The effect of 
the grid sensitivity, the time step, and an artificial dissipation have to be assessed along 
with the capability of existing turbulence models to capture unsteady flow physics and the 
transitional boundary layer. The Euler boundary layer procedure developed by Fan and 
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Lakshminarayana (1996) is very efficient due to the parabolic nature of viscous layers. 
But, this procedure is restricted to thin unseparated viscous layers, and its accuracy 
depends on the accuracy of the Euler solution. The Navier-Stokes code, on the other 
hand, is more general and does not depend on inviscid/viscid uncoupled procedure. But its 
disadvantage is in large CPU time due to large number of grids (typically y + =u T y/v=l) 
required for the computation of the amplitude and phase angle of various flow properties 
inside the viscous layers. 

1.1.3 Film Cooling 

The theoretical limit of the efficiency of a thermal power plant is governed by the 
Carnot cycle. The temperature of the heat sink is usually equal to or higher than the 
ambient temperature, thus the only way to improve the thermal efficiency is to increase the 
temperature of the heat source. An increase in the inlet turbine temperature is one of the 
most efficient means of advancing efficiency and weight characteristics of turbines. In gas 
turbines, the relation between the inlet turbine temperature and turbine efficiency is 
complex, and includes the compressor pressure ratio and bypass ratio as well as other 
parameters. 

Combustion of modem gas turbine fuels can provide stochastic temperature in 
excess of 2200 K, while modem materials cannot stand temperatures higher than 1200 - 
1400K. If the cycle temperature exceeds these values, blades should be cooled. The most 
common techniques for cooling are the convection cooling and the film cooling. For the 
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engines with turbine inlet temperature in excess of 1600K, only film cooling or a hybrid 
convective and film cooling techniques can provide acceptable surface temperature. 

A significant amount of research, predominately experimental, has been done 
during the last several decades to improve the efficiency of film cooling and to understand 
the aerothermal flow physics associated with the process. Early experimental research 
were carried out using simplified experimental conditions such as injection from a single 
hole in the direction of the flow on a flat plate. More advanced research utilized 
configurations of practical interest. The main problem associated with these efforts is in 
making accurate heat transfer measurements in a real machine. The influence of different 
characteristics, such as, blowing ratio, hole shape, injection angle, and turbulence intensity 
of the freestream on the film cooling effectiveness were investigated by Goldstein et al., 
(1987), Bergeles et al.,(1977), and Pietrzyk et al., (1990). Tekeishi et al., (1991) measured 
the film cooling effectiveness on the rotating turbine stage. Abhari and Epstein (1994) 
measured the time-resolved heat transfer on the rotor of a transonic turbine stage. A 
review of some works can be found in Margason (1993), Lakshminarayana (1996). 

Early attempts to predict film cooling effects were based on the parabolic or the 
partially-parabolic equations. Craward et al., (1980) used a boundary layer code to predict 
laterally averaged film cooling. Bergeles et al. (1977) used a semi-elliptic code with the 
prescribed constant velocity at the jet inlet. While codes based on partially parabolic 
equations and especially boundary layer codes are extremely effective, the predictions are 
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at best qualitative. Fully coupled procedures should be used to improve the accuracy of 
the prediction. 

Almost all current efforts in numerical modeling of the film cooling are based on 
coupled solutions. Current numerical efforts can be divided into two groups. The first 
group of researchers attempt to simulate simplified geometries such as a flat plate (Lylek 
and Zerkle (1993)) or a magnified model of the leading edge cooling (He et al. (1996) ). 
This approach provides good numerical resolution of the jet structure and is aimed at 
resolving the physics of cooling jet-mainstream interaction. Another approach is to 
simulate the flow in the real turbine geometry. Numerical simulations of film cooling flow 
in a turbine were performed by Hall et al.(1994), Vogel (1996), Garg and Abhari (1996). 
Due to the memory and CPU time limitation, only a limited number of grid points were 
distributed inside the film cooling hole. Lack of an adequate grid density reduces the 
accuracy of the prediction. 

The prediction of the film cooling is in its infant stage. Even though recent 
attempts are promising, none have been able to predict the film cooling effectiveness and 
jet-mainstream interaction accurately. This is mainly due to the numerical inaccuracy, 
turbulence model, and grid sensitivity. These issues will be addressed in this thesis. 

1.1.4 Three dimensional flow in turbine including tip leakage effect 

Most of axial turbomachines have a small clearance between the rotor blade tip 
and the casing. The presence of the tip gap generates the tip leakage flow, which has a 
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profound effect on the stage aerodynamics, efficiency, and noise vibration. According to 
Schaub et al. (1993), in a modem high performance high-pressure turbine, up to 30% of 
losses can be attributed to the presence of the tip leakage flow. Tip leakage flow results in 
reduced loading of the blade. Another significant effect is the modification of the heat 
transfer pattern due to the interaction between the tip leakage flow and the mainflow. The 
interaction between the tip leakage flow and the succeeding blade row results in an 
additional source of unsteadiness. Chopping of the tip vortex by the leading edge of the 
downstream blade produces turbulence and mixing, contributing to increased losses. 

The leakage flow has a complex three-dimensional structure. Development of the 
tip leakage flow is characterized by the complex interaction between the passage 
secondary flow, tip clearance vortex, blade wake, and the the endwall boundary layer. 
Significant efforts have been made to improve an understanding of the tip clearance and 
secondary flows in turbine. A comprehensive review of the experimental and 
computational research in this field can be found in Sjolander (1997). 

Due to the complexity of flow measurements, most of the research work was 
limited to the cascade flows; Langston et al. (1977), Gregory-Smith et al. (1988), and 
Yamamoto (1989). It is only recently that the emphasis is placed on the experimental work 
in real turbine stage configurations. Experimental measurements in actual rotors, (Joslyn 
and Dring, 1992 and Ristic et al., 1998) indicate significantly different leakage and 
secondary vortex structure in comparison with those observed in a linear cascade. One of 
the main conclusions of the experimental investigation by Ristic et al. (1998) is that the tip 
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clearance vortex is confined to the suction side corner of the blade, unlike in a cascade 
where the considerable pitchwise flow transport is observed. 

Several models based on the inviscid consideration of the tip leakage vortex (e.g., 
Lakshminarayana, 1970; Senoo and Ishida, 1987) are successfully used by the turbine 
industry in their design systems. However, further efficiency improvement requires a better 
understanding of the complex secondary vortex structure including the analysis of the loss 
origin. Utilization of the numerical modeling is a valuable tool in the achieving this 
objective. 

1.2 Objectives and thesis organization. 

Many aspects of the turbomachinery flow physics are still unresolved. Further 
progress can be achieved through a systematic application of the computational technique 
to the investigation of turbomachinery flows. Extensive validation and certification process 
is a necessary step in order to establish confidence in numerical simulations. The main 
objective of the thesis is to contribute to a better understanding of the turbomachine 
aerothermodynamics through the development and utilization of the numerical modeling, 
with special emphasis on the code validation and calibration, turbulence and transition 
modeling aspects. 

The main objective of the research is achieved through the accomplishment of the 


following tasks: 
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(1) Development of the unsteady Navier-Stokes solver, which can provide an 
efficient and flexible modeling of turbomachinery flows. This task includes an 
improvement in the computational efficiency of the solver for unsteady applications (dual 
time step approach) and improvements in the code flexibility and physical models (i.e. 
incorporation of the multiblock, flexible boundary conditions, incorporating of wide range 
of turbulence and transition models etc.). 

(2) Establishment of the solver’s reliability range for turbomachinery unsteady 

flows. 

(3) Application of the solver to the investigation of complex turbomachinery flows 
in order to gain a better understanding of flow physics. This task includes the analysis of: 

(3.1) Unsteady transitional boundary layer: Assessment of the turbulence 
models for their ability to simulate wake induced transition. Analysis of the effect of the 
unsteady transitional boundary layer development on the turbomachinery performance. 

(3.2) Different modes of steady transition varying from bypass transition in 
the attached flow to the transition over a laminar separation bubble. Effect of the 
numerical scheme on the accuracy of the prediction. 

(3.2) Rotor-stator interaction in the form of the upstream wake transport 
through the stage. Analysis of the upstream wake decay mechanism (e.g. contribution of 
inviscid stretching and viscous dissipation), and its effect on the unsteady flow loss 
generation. 

(3.3) Flow physics of the rotor-stator interaction at off-design conditions. 
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(3.4) Analysis of secondary flows due to the jet-mainstream interaction in 
film cooling configurations. Identification of the sources of the aerodynamic and heat 
losses due to the presence of the vortices. 

(3.5) Analysis of the secondary flow in a turbine rotor, including the tip 
vortex - passage vortex interaction. Effect of the tip clearance flow on the rotor efficiency. 
Development and decay of the tip clearance vortex. 

The main steps of the research presented in this thesis are illustrated in Fig. 1-1. 

1.3 Contribution of the thesis 

An unsteady compressible Navier-Stokes code based on the pseudo-time 
acceleration technique has been developed. Incorporation of the pseudo-time approach has 
enabled efficient unsteady simulation with CPU utilization improvement from 5 to 25 
times in comparison with the original code. An analysis of the scheme has been carried out 
to assess different approaches to the discretization of the time derivatives in the pseudo- 
time based scheme. Results of this analysis have been used to establish correction for the 
local pseudo-time step (iterative parameter) to provide efficient and stable unsteady 
calculations. Multiblock feature has been added to the code in order to simplify grid 
generation process and improve grid quality for cases with complex multidomain 
configurations such as the film cooling and the rotor-stator interaction. 

Extensive validation of the code has been performed to assess sensitivity of the 
solver to grid characteristics and artificial dissipation for the complex turbomachinery 
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flows (e.g., preservation of the accurate decay of the unsteady wake, the accurate 
development of the unsteady boundary layer etc.) 

The code has been used for the investigation of the unsteady flow physics in 
turbomachinery blade rows. Analysis of the wake induced transition on a turbine and a 
compressor blade has been carried out. A number of turbulence models have been 
assessed for their ability to provide an accurate prediction of the rotor-stator interaction 
effects, including the wake induced transition. 

Detailed simulation has been performed to investigate the transport of the 
upstream wake through the turbine rotor and the mechanism responsible for the wake 
decay. The contribution of different physical mechanisms; potential interaction viscous 
dissipation, and inviscid stretching has been analyzed. The viscous dissipation has been 
found to be a major contributor into the overall wake decay. However, the wake 
stretching inside the blade passage is predominantly inviscid. Distribution of losses 
correlates with this conclusion, as most of the losses due to the unsteady interaction are 
concentrated upstream of the leading edge (wake mixing losses). 

The flow in a low pressure turbine at different Reynolds numbers and freestream 
turbulence levels has been studied. Variation of the flow condition results in a different 
type of transition, varying from a bypass transition to a separated flow transition. Different 
approacheas to the transition modeling (different turbulence and transition models) have 
been assessed to establish their ability to predict the transitional flow in a low pressure 
turbine within the range of flow conditions encountered in practice. The analysis 
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performed shows that the prediction of the flow with the transition over a separation 
bubble is more sensitive to the numerical aspects of the scheme in comparison with the 
attached flow transition. 

Investigation of the complex aerothermal field due to the leading edge at a 
compound angle has been carried out. Results of the modeling have been used to analyze 
the vortex structure associated with the coolant jet-freestream interaction to understand 
the effect of different vortices on the cooling effectiveness and aerothermal losses. The 
system of vortices has been found to be essentially different from those observed in a flat 
plate configurations. Effect of the inlet turbulence and the length scale on the aerothermal 
field has been examined. It has been found that the inlet turbulence scale has a profound 
effect on the accuracy of the prediction. This influence is significantly higher in 
comparison with those observed in the boundary layer flow due to the intense mixing and 
entrainment of the ambient fluid into the coolant jet vortex structure. 

Numerical modeling is employed to gain a better understanding of the secondary 
flow in the Penn State rotor. Result of the simulation is used to interpret tip vortex 
development, its interaction with secondary flow and vortices. Secondary flow vortex 
structure is analyzed to estimate its contribution in overall loss generation. 

1.4 Thesis organization 

The governing equation and numerical scheme employed are described in Chapter 
2. Code development includes two major components; the development of an efficient 
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unsteady solver based on a pseudo-time step approach and incorporation of the multiblock 
capabilities. The results of the stability analysis of the dual time step scheme are used to 
establish the optimum pseudo-time step based on local CFL, Von Neuman numbers and 
the ratio of pseudo-to-physical time steps. 

The application of CFD analysis to the investigation of the physical problem 
requires the establishment of the reliability and the accuracy of the code. Results of the 
code validation and verification are presented in Chapter 3. The emphasis is on the 
assessment of flow features crucial for the numerical model developed in the thesis. The 
influence of the numerical aspects (artificial dissipation, grid density etc.) on the accuracy 
of the prediction of the ffeestream wake propagation, transition, flow over a cylinder is 
analyzed. Based on the results of the test cases, modifications to the k-e model to improve 
the accuracy in regions with dominant normal stresses are discussed. 

Investigation of the unsteady flow in a compressor stage is presented in Chapter 4. 
This research is carried out in three major aspects. First aspect is the validation of the code 
against the experimental data and the establishment of the number of pseudo time and 
physical time steps required for the accurate simulation of the unsteady flow. Another 
aspect considered is a comparison of the current prediction with the prediction based on 
the Euler/boundary layer approach from the point of view of accuracy and efficiency. 
Development of the unsteady boundary layer, including the unsteady transition zone, as 
well as the upstream wake-profile wake interaction effects is also discussed. 
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Chapter 5 is aimed at an improved understanding of the flow physics in turbines 
through the integrated computational and experimental study. The prediction has been 
validated not only against the blade surface experimental data (i.e., unsteady surface 
pressure distribution), but also against instantaneous blade-to-blade velocity acquired from 
an LDV. The sources of additional losses due to the unsteady interaction are analyzed. 
The physics of the upstream wake transport and decay are investigated. The emphasis is 
on the contribution of different mechanisms responsible for the overall wake decay. The 
results presented on the wake induced transition show the ability of the code to simulate 
major features associated with the unsteady transition, with the exception of the calmed 
region. 

Chapter 6 incorporates research efforts on the numerical simulation of the flow in a 
low pressure turbine. A range of parameters are considered. These variations correspond 
to different types of transition, from bypass transition in the attached flow to the transition 
over a laminar separation bubble. Assessment of different turbulence and transition models 
is presented. Predictions of the transitional flow based on k-e and ARS turbulence models 
are compared with the prediction based on the utilization of the transition model. 
Extensive evaluation of the effect of the artificial dissipation on the accuracy of the 
transition prediction is carried out to estimate the accuracy range. 

Jet-main flow interaction may lead to the generation of a vortex structure and, 
consequently, to additional aerodynamic losses and heat transfer. Examples of this type of 
flows are the leading edge film cooling considered in Chapter 7, and the tip clearance flow 
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in turbine, discussed in Chapter 8. The emphasis of this research is on an achievement of a 
better understanding of the vortex structure due to the jet-mainflow interaction, including 
the identification of sources of vortices and associated losses (both aerodynamic and 
thermal). Comparison of the simulated flow of the leading edge film cooling model with 
the experimental data shows the ability of the code to predict complex vortex structure 
associated with the film jet-main flow interaction. Results of the simulation are used to 
gain better understanding of factors affecting thermal and aerodynamic efficiency of the 
leading edge film cooling. Numerical analysis of the tip leakage flow in a turbine is utilized 
to investigate secondary flow physics in the rotor, including interaction between secondary 
vortices, tip leakage vortex, and rotor wake. 

Conclusions from the current research as well as suggestions for the future 
research are summarized in Chapter 9. 
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Chapter 2 


GOVERNING EQUATION AND NUMERICAL PROCEDURE 

A three-dimensional steady/unsteady Navier-Stokes solver has been developed in 
this thesis. Present development is based on the extension of the original solver developed 
earlier (Kunz and Lakshminarayana, 1992). Pseudo-time stepping has been incorporated 
to enable an efficient unsteady computation. Multiblock feature has been incorporated to 
make the code more flexible for the computation of the flow in complex topologies. A 
description of the code development as well as the numerical procedure, and turbulence 
models utilized are described in this chapter. 


2.1 Governing equations and numerical procedure 


Applying the Favre averaging procedure to the continuity, momentum, and energy 
equations, the five mean flow equations can be written in Cartesian tensor form as : 

^p+3p-(p(7 ) = 0 
axj 

^^^ + ^(pU i U J ) = -^ + ^-r i j-p£ i j k co j £ ldm co l x m -2p£ ljk co J U k 
dt OX : OX; OX; 




a(P ‘°' ,) + + ^Wj) = #-(t7 ,t, - q,) 


dt 


dx, 


Where: / = - - f fdt - ti 


time averaging of / 


[ 2 - 1 ] 
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/ = — - Favre averaging 
P 

T y = T Uj — pu*u' - effective stress tensor 
<Ji — q ti + P u i e - effective heat flux vector 

e 0R = p (e + W 2 12- (O 2 r 2 12 ) -energy transport variable, assuming rotation 

vector is coincident with x-axis and r is a distance to the axis 

Reynolds stresses and heat flux components are calculated using the eddy viscosity 

hypothesis or higher order turbulence closure. 

For the stability analysis of the numerical scheme presented below equations [2-1] 

can be rewritten in a matrix form: 

dQ d(E'+E v ) d(F,+F v ) d(G i +G v ) 

= 1 1 ho 

dt ox oy oz 

12 - 2 ] 

e=(p pc. pc, pc. P e oR ) ~ vector of conservative variables 
E t , F t ,G t - inviscid flux vectors, 

E V ,F V ,G V - vectors of viscous terms, 

S - source vector. 

Explicit four-stage Runge-Kutta scheme is used for time integration of both main- 
flow and turbulence equations. A compact second order accurate central difference flux 
evaluation scheme is employed for the convection terms. Diffusion terms are discretized 
using second order accurate central differences. For the mean flow equations a fourth 
order artificial dissipation is included to damp high wave number errors. Second order 
dissipation is used to improve the shock capturing. Eigenvalue and velocity scaling are 
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used to optimize the amount of the artificial dissipation. Multigrid and an implicit residual 
smoothing are used to improve the convergence characteristics of the steady solver. 

Numerical simulation of the unsteady flow requires special efforts to reduce 
possible reflection at the boundaries. One and two-dimensional, Giles’ type, non-reflecting 
boundary conditions are incorporated to minimize the reflection at boundaries and to 
minimize the computational domain (Fan and Lakshminarayana, 1996). 

2.1.1 Turbulence closure 

Turbulence equations are descretized in a manner similar to those for mean flow 
equations. "Lagged” approach is utilized for the computation of turbulence equations, i.e., 
k and £ values at previous time step are used to calculate the eddy viscosity at the current 
step. The presence of the source term in the turbulence equation results in a stability 
problem during the initial convergence period. Two mechanisms are used to ensure a 
stable calculation; utilization of the underelaxation factor, <J>, for k and £ equations in 
addition to the time step based on the mean equation and enforcement of an eddy viscosity 
limit. The maximum ratio of iVl-ii is set equal to 10-100 during the initial convergence 
period with yhe further increase to 1000-10000 to ensure the correct solution. It was 
found that the utilization of <{>=0.6 in the case of two-dimensional flow and 4>=0.75 for the 
three-dimensional flow improves convergence characteristics of the solver. 

2.1.2 Two-eauation models 

In eddy-viscosity models effective stress tensor and effective heat flux vector are 
defined: 
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T ij =r uj ~ P u ' u * and q t - q u + pu t e 


Reynolds stress is calculated from 


- P«,« y = P, 


dx , dx i 3 dx k 


-f a « p; 


heat flux component is given 

„ = _ r ElK 

* ' Pr, a*, 


[2-3] 


Turbulent eddy viscosity is computed using Prandtl-Kolmogorov relation 

c f pk 2 
1 £ 

[2-4] 

Transport equations for turbulence kinetic energy and turbulence dissipation ratio (for the 
simplification, in subsequent equations averaging symbols are omitted) 


jw t djpUjk) _ a 

dt dx, dx , 


J L 


, P, x M 
(M,+—h 


<7* dxj 


+ pp k - pie + D) 


dips) | dipUj£) _ a 

9/ 9jc. dx 


> L 


/ P, \ de 
(p,+— ) 


tr £ 9x yJ 


+ P , fc'tl/lP* PP 

k 


[2-5] 


[2-6] 


Here f t , f 2 , fj,, c^, c c i,D, E are low Reynolds number functions and constants 
described below. 


n " * 

P k = —UjUj 


W L 

dXj 


- production of turbulent kinetic energy 


A number of low-Reynolds number k-e turbulence models; Chien (1982), denoted 
as CH, Lam-Bremhost (1981), denoted as LB, and Fan-Lakshminarayana-Bamett (1993), 
denoted as FLB, are utilized for the turbulence closure. A summary of the constant and 
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near wall function for different turbulence models is given in Table 2-1 and Table 2-2. The 
solution of the turbulence transport equations is numerically coupled with the solution of 
the main flow equations. A second and fourth order artificial dissipation is included in the 
turbulence transport equations. 


Table 2-1 Low-Reynolds-number functions used in turbulence model 


Model 

Code 


Chien 

Ch 

l-exp(-0.01 lSy*) 

Lam-Bremhost 

LB 

[ 1 -exp(-0.0 1 65Rey)] z ( 1 +20.5/Ret) 

Fan-Lakshminarayana- 

Bamett 

FLB 

0.4fJ^/Re, +(\-0AfJ ^]Re t )[l-exp(-Rey/42.63)] 3 


Table 2-2 Low-Reynolds-number functions used in turbulence models 


Code 

f. 

f 2 

D 

E 

CH 

1.0 

l-0.22exp(-Re t 2 /36)) 


-2v(e/y 2 )exp(-0.5y + ) 

LB 


l-exp(-Ret 2 ) 

0 

0 

FLB 

1.0 

1 -2/9exp(-Re t 2 /36))f w 2 

0 

0 


where the near-wall function in the FLB model is given by 


L = l-exp 




2.30 


‘.(V s 


Re. 


2.30 8.89 


Re v 

1 - exp( 

K 20 


-13 
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Constants Q=0.09,Cei=1.44,Ce2=1.92,a k =1.0,a E =1.3 in the LB and FLB models 
are the same as those used in the standard (high-ReynoIds-number) k-e model given by 
Launder and Spalding (1974). Chien's k-e model has slightly different values for 
C e ,= 1.35,0^=1.80. 


2.1.3 Algebraic Reynolds Stress Model 

Numerical simulation based on the first order turbulence closure can be applied to 
a wide range of cases and may improve the accuracy of the prediction in comparison with 
algebraic turbulence models. However, deficiency of these models associated with the 
Boussinesq approximation and empirical correlation used to derive these models leads to 
less precise solution in the case of complex flows (Lakshminarayana, 1986). Non- 
equilibrium flows, flows with streamwise curvature and rotation, flows with injection (e.g., 
film cooling) are examples when the first order turbulence closure does not provide an 
adequate level of accuracy. Flow computation based on a second degree closure and a 
subgrid turbulence modeling (LES) demonstrated the potential for the improvement in the 
turbulence flow prediction. Complex models generally require more CPU time. Another 
factor affecting the wide acceptance of more complex models is a potentially increased 
dependence on the numerical stability. More complex structure may lead to a less robust 
and, as a result, less reliable prediction. It is more difficult to develop a stable code in the 
case of full Reynolds Stress (FRSM) models. To overcome stability limitations many 
FRSM solvers utilized a simplified approach for the near wall region (wall function, one- 
equation models), thus decreasing the accuracy of the flow resolution near the wall. 
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Rodi, 1976, suggested a simplified algebraic expression for the component of the 
Reynolds stress tensor. This model is based on the assumption that the transport of 
Reynolds stress components is locally proportional to the transport of the turbulent kinetic 
energy. ARSM may be considered as a compromise between the two-equation and higher 
order models. Implementation of the ARSM does not lead to a significant increase in the 
CPU time and requires an inversion of the 6X6 (multidimensional case, implicit ARSM). 
ARSM uses the following expression to calculate Reynolds Stress component: 


- u'u' = -k[R tJ (2 - C 2 ) / 2 + (^ - 3)0 - C 2 )]/ ^P kk + £(C, - 1) 




where: =P k = -u'u' - production of k 

d Xj 

P r = -u'u' — — i- + u'u' - production of Reynolds stresses 

,J ‘ dx k 1 dx k 

R, = -2 03 p (£ ipk u'u k +£ Jpk u'u' k ) 

Ci=1.5, C2=0.6 - model constants 


--8,jk 

3 ' 


[2-7] 


In current research, a hybrid model is utilized in a near wall region. Laminar 


sublayer and overlap region are calculated using k-£ equation. Matching function based on 
Re y is incorporated to smooth the transition between regions calculated using k-e and 


ARS models: 


Re 


/* 2 (- 


tanh (B — —I y* match - 1) 
H 1.831 


tanh(j8) 


+ 1 ) 


where: (3 is a slope constant 
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y + match - matching location 

Kis = fm R. ARSM 0 — f m )^k-e 

2.2 Development of the multiblock solver 

The quality of the grid used in the numerical simulation has a serious impact on the 
accuracy of the numerical simulation. In the case of the complex computational topology, 
it is appropriate to divide the computational domain into sub-domains and consider these 
sub-domains as separate computational blocks. Utilization of the multidomain structure 
enables better optimization of the grid point -distribution, thus leading to an improved grid 
quality. The multiblock approach is also useful in single-connected regions with relatively 
simple topology, such as a turbomachinery blade row. In this case, different grid types can 
be more suitable for different regions. A C-type grid is better suited for the near blade 
region including the leading edge, while an H-type grid can be utilized for a mid passage 
and inlet/outlet regions. The multiblock approach can also simplify the implementation of 
the zonal approach, i.e., application of different mathematical models in different regions. 
Numerical simulation of configurations, with relative motion is another case where a 
multiblock approach is very helpful. An example of this type of problem is the stator/rotor 
interaction. Simultaneous calculation of the flow field in several rows improves the 
accuracy of the numerical prediction. A multiblock approach, with grids stationary in 
relation to the corresponding blade row and moving relatively to each other, is the best 
approach for the computation of flows with rotor/stator interaction. 
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2.2. 1 Structure of the multiblock solver 

The three-dimensional viscous code has been extended to include a multiblock 

feature to enable a computation of the rotor/stator interaction and other complex 
turbomachinery flows. A number of alternative approaches to the development of the 
multiblock solver have been considered. The first approach is the extension of the existing 
array structure with an additional index representing a block number. This variant has two 
essential drawbacks. First, it leads to an excessive use of memory. The total required 
memory is equal to the memory required for the storage of the largest block times the 
number of blocks. Therefore, efficient memory allocation can not be achieved in the case 
of uneven block sizes. Another significant drawback is the need to rewrite the whole 
program, resulting in extensive additional debugging and incompatibility with previous 
version of the solver. Utilization of Fortran 90 array type v with variable element length has 
been the second considered approach. Preliminary tests have indicated that this may lead 
to a certain level of the performance degradation (about 15-20%). 

The replacement of multi-dimensional array structure with one-dimensional arrays 
enables an efficient memory allocation. However, it still requires a significant rewriting of 
the code. A multilevel approach has been adopted in the current research to combine 
advantages of one-dimensional arrays with the original code preservation. The code is 
divided into three modules: “control”, “communication” and “kernel” (Fig. 2-1). The main 
purpose of the control module is to provide the switching between “kernel” and 
“communication” modules: 


DIMENSION X (n«ieJ , . 
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DO block_number =1 , number_of_block 

CALL kernel (id (block_number) , jd (block_number) , 
kd (block_number ) , X ( i_lock (block_number) , . . . ) 

CALL communication (X, . . . . ) 

enddo 

SUBROUTINE kernel (i di m/j dim/ kdim/X/ • . • ) 

DIMENS ION X ( idim/ j dim/ k dim ) 

end ! kernel 

SUBROUTINE communication (X, . . . ) 

DIMENSION X ( n e iem) 

end (communication 

In “control” module all data is represented as one-dimensional array structure. 
“Kernel” module interface works as a switch between one-dimensional and 
multidimensional data representation. Fortran passes an actual argument of subroutine as a 
reference. i_lock (block_number) corresponds to the first element of array X, 
belonging to block_number. The specification of X (...) as actual argument is equal to 
the employment of a pointer to indicate the position of corresponding block in X array. 
Inside the “kernel” module there is no information about existence of other blocks. Same 
routines as in the original single block version of the code are used to perform 
calculations. “Communication” module is the only module with the simultaneous access to 
the elements belonging to different blocks. Interblock data transfer requires knowledge of 
the position of the interface elements in one-dimensional array structure. This task is 
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performed during the preprocessor stage based on a given grid and interblock topology 
description. 

The current structure of the multiblock is well suited for the parallel computer 
implementation. The development of the Message Passing Interface, MPI, based version 
of the code for parallel computers requires the modification of three locations only. 
“Kernel” routine calls should be set to be executed in parallel. In “communication” 
module, MPI send/receive routine should be added to provide the data exchange between 
elements stored in memory blocks situated at different CPU. Utilization of pointers 
makes it easy to move the code from distributed memory systems to shared memory 
systems. Based on the system type, the preprocessor can calculate i_lock suitable for 
the corresponding system. For example, for the distributed memory system n e i em can be 
set equal to the number of elements in the largest block, while 
i_lock (block_number) will be equal to one for all blocks. For the shared memory 
system n elem is equal to total number of elements in all blocks. 

The information exchange between blocks is a crucial component for the 
successful development of the multiblock solver. The communication procedure must be 
efficient, robust and be able to preserve conservation properties. Overlaid and patched 
grids are the most common types of multiblock grids. The current solver uses overlaid 
grids. At each interface grids are overlapped by one grid point. Data assigned to this point 
is based on the solution in the adjustment block. To achieve an accurate preservation of 
conservative properties a conformal interface has been chosen as a primary mechanism for 
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the information transfer. Conformal interface requires that the boundary ghost cell must be 
coincident with the inner cell of the adjustment block. No interpolation is required. Data is 
directly passed from one block to another. For some computational topologies conformal 
interface may complicate grid generation. To increase a flexibility of the code, the second 
type of the interface with data interpolation has been added. However, all multiblock 
results presented in the current thesis are based on the conformal interface approach. 

The explicit nature of the numerical solver utilized in this research limits the 
distance information can propagate during each iteration. Thus, the separation of the 
computational domain into subdomains does not affect convergence characteristics. 
Interblock interfaces (in the case of the conformal interfaces) are essentially invisible for 
the solver. Simulations based on multiblock and the single block configuration (assuming 
that the multiblock topology can be represented as single block) result in identical 
solutions. Convergence behavior is also practically identical. The only exemption is the 
flow with rapidly changing conditions across the interblock interface. This is due to the 
fact that an artificial dissipation is calculated using one-side finite differences at each side 
of the interblock interface. One side differences are used because one cell overlap does not 
provide enough grid point to calculate the 4 th artificial dissipation term. 

In comparison with the original version, boundary conditions are treated point-by- 
point in the multiblock version. This enables simulation of any general configuration (Fig. 
2-2). Point-by-point boundary conditions lead to an extensive calling of small subroutines 
and the problem with vectorization of this part of the code. The first problem can be 
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overcome by the inlining during the compilation. According to test cases CPU, overhead 
due to non-vectorization is less then 1-2%. 

Developed multiblock version of the solver has been utilized for the computation 
of the multidomain film cooling configuration presented in Chapter 7. Implementation of 
the multiblock approach has simplified the grid generation process and has improved the 
quality of the computational grid. 


2.3 Pseudo-time Acceleration 


Time-marching schemes are one of the most widely used methods for the 
numerical simulation of the compressible Navier-Stokes equations. Even though these 
schemes approximate unsteady equations, discretization errors and acceleration techniques 
can totally destroy time accuracy. During a steady code development the main emphasis is 


done on the minimization of the amplification matrix amplitude 



in order to 


increase the convergence rate. Phase error does not play any significant role. In contrast, 
time accurate computations require that the amplitude of G should be close to unity for all 
harmonics to be resolved. Phase angle error also has a major influence on the time 
accuracy. Thus special efforts should be undertaken in order to apply time marching 
schemes, developed for the steady state calculation, to the unsteady numerical simulation. 
One of the additional limitations is the restriction on a time step, i.e., the time step should 
be constant for all cells to preserve temporal accuracy. When an explicit, time marching 
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code is employed for the unsteady computation, the time step is limited by the minimum 
time based on stability considerations. Due to this, the actual time step is much smaller 
than the time step needed to achieve required temporal accuracy, especially in the case of 
highly stretched grids in viscous flows. Furthermore, it is impossible to use common 
acceleration techniques, such as multigrid and implicit residual smoothing. These 
techniques generally affect the temporal accuracy of computations. These limitations lead 
to a large increase in CPU time utilization for unsteady viscous computations on highly 
stretched grids. 

Implicit schemes do not imply direct limitations on a time step. However, most 
implicit schemes introduce additional linearization error, especially if special technique 
(e.g. approximate factorization) is used to simplify matrix inversion. The amplitude of the 
error is proportional to the utilized time step. Thus, implicit schemes developed for steady 
solvers do not provide the efficient and accurate simulation for unsteady problems also. 
These difficulties can be overcome through the introduction of dual step calculations. 

For example, the governing equation [2-2] for two dimensional time marching 
problems can be written in the form: 



Where R(Q) = - 


d(E,+E v ) d(F l+ F v ) 
dx dy 


[ 2 - 8 ] 


Original scheme uses one iteration to obtain the solution at a new time level, t n+l = 
t n +At. Iteration parameter t plays a role of the physical time step. In a dual step or pseudo- 
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time step approach physical time derivative ^^-is considered as an additional term in the 

at 

equation. To obtain the solution at t n+l = t n +At iterations are performed in mathematical 
space x rather than in physical temporal space: 

Equation [2-8] can be rewritten in the form: 


dQ _dQ 
dx dt 


+R(Q ) 


or d S- = R\Q) 


dx 


[2-9] 

Physical derivative — can be discretized in an implicit manner, thus removing 
dt 


limitation on the physical time step. There is no direct influence of the solver 
characteristics on the temporal accuracy. The original steady solver may be used as an 
iterative procedure during the inner cycle to achieve the time accurate solution at t n+1 = 
t n +At. The only effect of the discretization scheme (temporal in a pseudo time space) on 
the temporal accuracy is the level of convergence during the inner cycle. 

Equations [2-8] and [2-9] are similar, the only difference is the presence of 
additional terms in the residual. Steady computational codes which were developed to 
solve the equation [2-8] can be applied to the equation [2-9] with minor modifications. 
The same acceleration techniques that are used in steady state calculations can be utilized 
in the inner cycle. However, the presence of the source term and the rate of convergence 
during inner iterations can affect the efficiency of the code and should be carefully 
analyzed. The form of discretization used for physical temporal derivatives is the factor 
affecting the stability and convergence rate of the internal cycle. This factor should be also 


taken into consideration. 
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Implicit schemes with approximate-Newton methods can be considered as a 
subgroup of schemes with a dual time step. Different schemes with the dual time stepping 
have been successfully used to simulate unsteady problems during recent years (Amone et 
al. (1995), Hall (1995), Alonso etal. (1994) Daily etal. (1995)). 


2.3.1 Stability analysis of the scheme with pseudo time acceleration 

Different approaches have been used in the discretization of the physical derivative 
in pseudo-time step techniques. These approaches for the discretization of the equation [2- 
9], can be expressed as follows: 

(l+a' k ^c)A Q k =-a k At(S k ~ l +R k ~ l )+a' k ^cAQ* -1 

[ 2 - 10 ] 

where: a k - coefficients of Runge-Kutta scheme 

a[ al - auxiliary coefficient of Runge-Kutta scheme with pseudo time stepping 
S - additional source term due to the pseudo-time stepping 

Second-order accurate discretization of time derivatives can be written 


30 3 o k ~ l —4 o n +0 n_l 

as: — = S = — — — — , in this case constant c=2/3. Values of the coefficients 


dt 


2At 


used by various authors are given below: 

1) a' k =0,al = 0 ; explicit treatment of physical derivatives in internal cycle 


(Amone et al. (1995) and Hall(1995)), denoted as scheme 1. 

2) Q^=l,a' = 0 ; implicit discretization of physical derivatives (Weiss et 


al.(1995)), denoted as scheme 2. 
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3) o£ =OC k ,c( =a k ; also implicit discretization of physical derivatives (Melson et al. 
(1994)), denoted as scheme 3. 

4) Oc[ =<X k ,a' k =0, denoted as scheme 4 

The presence of the additional source term in equation [2-10] changes the 
behavior of the scheme. Melson et al.( 1 994) carried out VonNeuman stability analysis of a 
five-stage Runge-Kutta scheme. This analysis has been applied to the four-stage scheme 
and has been extended to include different cases described above. The following two- 
dimensional model equation has been utilized for this analysis: 


du du du , . du ,d 2 u d 2 u . 3 d 4 u . 3 d 4 u 

dr dt dx dy dx 2 dy 2 dx 4 dy 4 

[ 2 - 11 ] 

where: a,b - model transport velocities 

Equation [2-11] is discretized using a four-stage Runge-Kutta scheme with 

following coefficients; ai=l/4 a 2 =l/3 a 3 =l/2 ou=l. A numerical scheme with the 

evaluation of all terms at all stages (denoted as scheme 1, scheme 2, etc.) as well as 

schemes with the evaluation of source and viscous terms only at the first stage (denoted as 

scheme la, scheme 2a, etc.) has been considered. An amplification factor of 

, aAr . bAr xAr tAr , LaAr , kJbAr s At , , . , 

has been denved 


using the symbolic computation program Mathematica. Only one-dimensional results are 
presented here for the sake of clarity. The second dimension does not principally alter the 
results, but makes the evaluation and interpretation more difficult. 
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The ratio of the pseudo-time step to the physical time step, 8 =At/Ai, plays a crucial 
role in the behavior of the scheme. The results of the stability analysis are summarized in 
Fig. 2-3-Fig. 2-8. An average amplification factor g 0 „ n is plotted as a function of the 

CFL number, VonNeuman number, and 8. This approach has been chosen instead of a 
more standard plot (real and imaginary parts of the amplification factor) to clarify the 
influence of various parameters. The plot is bounded by the stability surface with I g I = 1 . 
Explicit discretization of the physical derivative leads to a linear decrease in the maximum 
allowable VonNeuman number with the increasing 8 (scheme la. Fig. 2-3). On the 
contrary, an implicit evaluation of the physical time derivative (scheme 2a) results in an 
extended stability region at higher 8 (Fig. 2-4). The extended stability region indicates an 
advantage of the implicit evaluation of the physical temporal derivative. However, the 
distribution of g 0)r/2 suggests that within the stability region, a scheme with an explicit 
evaluation of the physical time derivative may possess better convergence characteristic as 
a result of lower g 0)r/2 . For both schemes, la and 2a, “optimal” X and a (at each value 

of 8) provide identical amount of error damping. Numerical modeling confirmed that 
scheme 2a and scheme la achieve similar convergence rate during the inner cycle if the 
“operational curve” (Fig. 2-3, Fig. 2-4) is used to adjust the calculation of the pseudo-time 
step t. These curves are used to adjust the time step according to the local value of 8. 

The presence of the additional source term results in another positive feature. The 
higher the ratio of the pseudo time step to the physical time step the more intensive 
decrease of the amplification factor is observed at low wave numbers. This is especially 


43 


beneficial for the numerical simulation of the wake propagation in a turbomachinery stage. 
At the beginning of each inner cycle an error spectrum is close to the spectrum of the local 
unsteadiness. Local unsteadiness is dominant by the frequencies associated with a relative 
rotor-stator movement. Usually only first few harmonics, based on the rotation passing 
frequency, have significant amplitude. These harmonics correspond to low wave numbers. 
Therefore, the utilization of the pseudo time stepping provides a very efficient mechanism 
of the error elimination in the region outside of the boundary layer (i.e., the zone with high 
5). Damping at low wave numbers is more profound in scheme la. Flat distribution of g at 
low wave numbers can be also used to explain the reduced efficiency of the multigrid 
acceleration when it is used in conjunction with a pseudo-time stepping. For the standard 
Runge-Kutta scheme the amplification factor is rapidly decreasing for low frequency 
errors. Calculations on a coarse grid effectively double a local wave number improving 
convergence. Practically constant distribution of the amplification factor at low wave 
numbers diminishes the effect of the multigrid acceleration in the case of a pseudo-time 
scheme. 

Figures Fig. 2-5 and Fig. 2-6 show an average amplification factor for schemes 1 
and 2 respectively. Evaluation of viscous and source terms at all stages extends the 
stability limit for both schemes. Low wave damping is also improved in scheme 1. 
However, simular to the steady state case, the advantage of the evolution of viscous and 
source terms at all stages is not significant enough to justify an additional CPU time 


required for this modification. 
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Last two schemes with the implicit evaluation of the physical time derivatives 
behave similar to those of scheme 2 (Fig. 2-7, Fig. 2-8). Scheme 3 possess more rapid 
decay of the g 0 „ n in comparison with the scheme 2. It is possible to expect that this may 

provide better convergence of the inner cycle for regions with 6-0.3- 1. An interesting 
feature of the scheme 4 is the independence of its stability limit from the ratio of pseudo 
time to physical steps. However, similar to all other schemes it provides better low 
frequency error elimination with increased 5. 


2.3.2 Artificial dissipation term adjustment for the solver with pseudo-time 
stepping 

Explicit treatment of the physical derivative in pseudo-time stepping imposes a 
limitation on the pseudo-time step during the inner iteration for the grid cell located in the 
middle of the blade passage. Thus the Courant-Fredrichs-Levy number for these grids may 
be significantly smaller than the maximum CFL^. Meanwhile, the artificial dissipation 
term is based on the local maximum CFL number. A one-dimensional simplified form of 
the governing equation [2-2] can be written as: 

^- + A^-+D(Q) + S(Q) = 0 
d T dx 

[ 2 - 12 ] 

where the artificial dissipation term : 


D(Q) = 


k 4 *S m Q*p(A) 
CFL* Ax 


[2-13] 


here p(A) - spectral radius of A 
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For CFL<CFLmax , D(Q) may become large even with a small variation in Q. To avoid an 
excessive level of the artificial dissipation , D(Q) was modified: 

D’ (Q)=D(Q)CFL/CFL max . 

[2-14] 

Incorporation of the pseudo-time acceleration has allowed an efficient simulation 
of unsteady turbomachinery flows, presented in subsequent chapters. Code based on 
pseudo-time approach requires from 5 to 30 times less CPU time in comparison with the 
original code. Scheme with explicit evaluation of the physical time step (scheme la) has 
been used, because it provides the same level of convergence as the scheme with the 
implicit evaluation of the physical time step (scheme 2a). 
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Fig. 2-1 Block scheme of the multiblock solver 
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Fig. 2-2 Allowable multiblock topology (in computational space) 
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Fig. 2-3 Stability limits and amplification factor distribution for scheme la 
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Fig. 2-4 Stability limits and amplification factor distribution for scheme 2a 
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Fig. 2-5 Stability limits and amplification factor distribution for scheme 1 
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Fig. 2-6 Stability limits and amplification factor distribution for scheme 2 
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Fig. 2-7 Stability limits and amplification factor distribution for scheme 3 
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Fig. 2-8 Stability limits and amplification factor distribution for scheme 4 


Chapter 3 


CODE VALIDATION AND MODIFICATION 

Code validation and verification are essential parts of the development of the CFD 
solver. Validation process of the current solver can be divided into two stages. First stage 
of the validation includes test cases to verify that the development and modification of the 
code did not introduce additional errors. The objective of the second stage of the 
validation process is to simulate test cases with complex features essential for accurate 
flow simulations in turbomachinery. In addition to results presented in this chapter, code 
validation against the existing experimental data for flow configurations discussed in 
following chapters is carried out to achieve confidence in the code and numerical 
simulations. 

3.1 Verification against analytical solution 

Numerical simulation of the inviscid, irrotational flow over a cylinder with inlet 
Mach number equal to 0. 1 has been carried out to verify accuracy of the code against the 
analytical solution. Flow over a cylinder possess a number of flow features which can be 
found in turbomachinery cascades, i.e., decelerating flow along the stagnation line near the 
leading edge etc. The existence of the analytical (potential) solution makes possible a 
quantitative comparison to assess the accuracy of the code. Utilization of the finite 
difference approach to the solution of the partial differential equations results in the 
procedure that effectively solves the equation, which is different from the original 
differential equation. Effect of the discretization procedure can be analyzed using partial 
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differential approximation (PDA) (Shokin and Janenko, 1985) of the equations. Deviation 
from the form of the original PDE may result in numerical errors, even if the computation 
could be carried out with the unlimited precision. To ensure stable calculation, most of the 
numerical codes include artificial dissipation terms. The absence of explicit artificial terms 
does not constitute the absence of the artificial dissipation but rather its increased value in 
the partial differential approximation. The analysis of PDA provides significant amount of 
information on the characteristics of the scheme, including the effect of the explicit and the 
implicit numerical dissipation. However, from the practical point of view, the verification 
and the validation of the code are more appropriate, even taken into the consideration the 
limited and empirical nature of the validation process. 

Flow over the cylinder has been simulated using 151 x 41, “H” type grid. Grid 
density is similar to the grid density (inviscid simulations and grids for viscous simulation 
outside the boundary layer) typically used for the computational analysis presented in 
following chapters. The comparison between the analytical solution and computational 
results (Fig. 3-1) reveals very good agreement for the computations with the appropriate 
amount of the artificial dissipation. Decelerating flow upstream of the cylinder is not very 
sensitive to the level of the artificial dissipation. Only use of very high values of let results 
in a deviation of the numerical solution from the exact one. How prediction downstream 
of the cylinder is not as accurate as the prediction of the upstream region. For cases with 
fourth order artificial dissipation coefficient Lt >0.015 a small separation zone develops 
near the downstream stagnation point. This phenomenon results in a flat velocity from x/R 
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= 1 to x/R = 1.05 (Fig. 3-1). Further increase of the artificial dissipation provides fully 
attached flow near the trailing edge. However, this results in a velocity field which has 
noticeably asymmetric downstream and upstream zones. The decreased total pressure 
indicates the unphysical energy dissipation. Quantitative analysis of the numerical error is 
presented in Fig. 3-2. First norm, || =max, y () represents the local error, while the 


second norm 


n i -n j 


is an indication of an overall accuracy of the prediction (n,,nj 


are grid dimensions). The computational error was calculated separately for the zone with 
x/R > 0 and x/R < 0. Stagnation points are the locations of the maximum errors for all 
cases with lc» < 0.04. Within a range of small values of k 4 = [0.05 - 0.015], the numerical 
error is practically zero outside a small zone near stagnation points. The development of 
the separation zone downstream of the cylinder results in fl-^ jump at Iq = 0.15. 

High artificial dissipation k > 0.04 leads to monotonic increase in H and ||| 2 

throughout the flow. This is an indication that at this level of the artificial dissipation, the 
type of the flow (accelerating or decelerating) as well as the flow grid alignment is 
irrelevant. Computations based on a doubled grid result in 30% decrease in || ( and have 


no significant effect on |J 2 for k 4 e [0.005, 002]. Based on this analysis, Li < 0.02 may be 
established as a requirement for the accurate flow simulation. 
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3.2 Influence of the artificial dissipation on the wake propagation and decay 

The correct simulation of the wake decay is an essential component of the rotor- 
stator interaction analysis. Excessive wake dissipation would result in the improper 
unsteady blade loading and unsteady losses. Numerical analysis of the freestream decay of 
the moving wake has been used to establish the requirement on grid characteristics and the 
level of an artificial dissipation. 

Moving wake has been simulated using nonuniform time dependent inlet boundary 
condition. Inlet total velocity was prescribed as U=U o +f(y-V 0W t). Time average flow angle 
is set at 45°. Function f is either Gaussian distribution or sine wave with different reduced 
frequencies. This configuration imitates the wake propagation in an axial gap in a relative 
frame of reference (without the potential effect). 

If inviscid flow model is considered, then any wake decay is due to the numerical 
dissipation. Numerical analysis has been carried out to establish the criteria (grid density, 
artificial dissipation) required to obtain an accurate wake propagation. 

Results of the numerical simulation are summarized in Fig. 3-3. The influence of the 
artificial dissipation on a wake decay is shown for two typical values of the fourth-order 
artificial dissipation coefficient, k». Twenty to thirty grid points per each wave width is 
necessary for an accurate prediction of the wake propagation. 

This flow can be analytically solved using the linearization procedure (see Appendix 
A): 

MS) = A, exp(-(27rOJ) 2 (-J- + (* 4 ^^1(1 + -U)£) 

Re n Ma 
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where £ = y-V HW x 

This formula correlates well with results of the numerical simulation and can be 
used to establish grid and requirement depending on the spectrum of the incoming 
unsteadiness. 

The results of this analysis are used as a guideline in generating grids for two 
additional test cases. The numerical solver was utilized to simulate the wake downstream 
of the flat plate measured by Chevray et al. (1969) (Fig. 3-4). The prediction is in good 
agreement with the data. The comparison of the numerical prediction of the far wake 
decay with the correlation due to Reynolds et al. (1979) for a cascade wake is also found 
to be in good agreement. 

Another objective of the current analysis is to verify the wake-outlet boundary 
condition interaction. Even though non-reflection boundary conditions are used in the 
current solver, certain amount of the wake damping occurs near the outlet boundary. This 
affect is limited by two-three grid points upstream of the outlet boundary. It does not 
generate any reflection wave. Thus, there is no adverse affect on the unsteady solution of 
near the boundary region. 

3.3 Steady turbulent and transitional boundary layer 

Numerical simulation of the transitional flow on a flat plate has been carried out to 
assess the ability of the code to predict the inception and the length of the transition. The 
test case chosen for this validation is T3a described by Savill (1992). The predictions from 
all three turbulence models (CH, LB, FLB) are compared with the data. A number of 
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investigators have used boundary layer codes to analyze the ability of the low-Re number 
turbulence models to predict the transition. This prediction is feasible since low-Re k-e 
equations model the transitional behavior of the boundary layer through low-Re functions. 
However, this prediction may not be very accurate, since low-Re functions employed are 
based on the fully turbulent flow. 

Uniform inlet flow and turbulence distribution are prescribed at the inlet, upstream 
of the leading edge. Numerical experiments indicate that the grid should be extensively 
stretched near the leading edge of the flat plate to minimize an effect of the singularity 
point and to ensure an accurate prediction of the transitional boundary layer. The skin 
friction coefficient distribution, shown in Fig. 3-5, indicates that the LB model shows the 
best agreement with the data for the low turbulence intensity, while the CH model predicts a 
very premature transition and the longest transitional length. Several ’numerical’ factors 
(artificial dissipation, grid density etc.) are found to have an appreciable effect on the 
prediction of the transitional region in comparison with laminar and fully turbulent zones. 
Possible variation of the Cf coefficient due to the variation of ’numerical’ factors is shown 
for the LB and FLB models. Numerical simulation of the flat plate flow with the higher 
inlet turbulence intensity, more typical for turbomachinery applications, shows that the LB 
model gives an earlier transition; while the FLB model is in better agreement with the data. 
This can be explained by the fact that the LB model is numerically, less stable than the FLB 
model. The momentum Reynolds number at the start and at the end of the transition, 
predicted from the Navier-Stokes code, is compared with the correlation of Abu-Ghannam 
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and Shaw (1980) and the boundary layer prediction (Fan & Lakshminarayana (1996)) in 
Fig. 3-6. The time-marching code predicts an earlier inception of the transition at low 
values of the turbulent intensity. Overall, the transition prediction by the Navier-Stokes 
code is close to the prediction by the boundary layer code based on the same turbulence 
model. 

In addition to the quality of the particular low-Re turbulence model, the level of 
the artificial dissipation plays a crucial role in the prediction of the transition region. For 
the configuration more complex than a flat plate flow, the potential error due to the 
presence of the artificial dissipation is even more serious. For example, for the transitional 
flow in Penn State turbine, the location of the transition and the skin friction coefficient 
beyond the transition, strongly depend on the value of Iq (Fig. 3-7). An excessive level of 
the artificial dissipation delays transition inception and ultimetely leads to the fully laminar 
boundary layer. An increase in ku results in the transition onset shifting from x/C*=0.65 to 
the trailing edge on the suction surface. Grid refinement reduced this dependency, 
however, as shown in Chapter 6, it is practically impossible totally eliminate the effect of 
artificial dissipation on the transitional prediction. The essential feature of the current 
solver is that for the calculations with L,< 0.015, the transition prediction is found to be 
independent of the level of the artificial dissipation coefficient. This value of the Ic* 
coefficient is found to be universal for other configurations ( the flat plate, compressor 
cascade, LP turbine), under the condition that the solution is otherwise grid independent. 
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3.4 Modification of the k-e model for the flows with dominant normal stresses 

Models based on the eddy viscosity do not provide an accurate solution for the 
flow that deviates from linear stress-strain relation. The problem is aggravated if two- 
equation turbulence model is applied to flows that are essentially different from those used 
for the derivation of the model coefficients. Flow with dominant normal stresses is the 
case where k-e model fails to provide a realistic prediction of the turbulence field. K-e 
model developed for the shear flow tends to overpredict the local level of turbulent kinetic 
energy produced by the normal stresses. A significant amount of work both experimental 
and computational were carried out to investigate these type of flow (e.g., Cooper et al., 
1993). Most of these efforts were concentrated on the analysis of the heat transfer 
associated with the impinging jet. Less attention to this problem was paid in aerodynamic 
simulations without heat transfer focus. This lack of attention can be explained by the 
relatively small influence of this problem on the overall accuracy of the prediction in many 
cases. Error in the energy redistribution between the mean flow and turbulence is about 
1 % for flows with Ma=0.4 and T u =8%. However, an excessive level of turbulence 
prediction seriously affects the transition development in a turbomachinery stage, even in 
the case of the utilization of the transition model. In the turbomachinery stage an 
overprediction of turbulence occurs at two major locations. First zone is the stagnation 
flow near the leading edge. The correct prediction of the turbulence kinetic energy at this 
location is especially needed if the boundary layer has the transition inception close to the 
leading edge (compressor cascade. Chapter 4). Adequate turbulence intensity near the 
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leading edge also makes the solver less “numerically” dependent (LP turbine. Chapter 6). 
The second region with the dominant normal stresses is the accelerated flow near the 
suction surface of the turbine blade. If no attention is paid to this region, turbulence 
intensity may be overpredicted by 3-7%, resulting in an earlier transition inception (turbine 
flows. Chapters 5 and 6). 

Application of more complex, in comparison with standard k-e, models (e.g., k-e- 
v 2 , Behnia et al., 1996) and second-moment closures (Craft et al„ 1993) showed an 
improved prediction of flows with dominant normal stresses. However, current analysis is 
limited to the modifications for a k-e model. 

3.4.1 Modification of the turbulence model for leading edge flow 

The modification of the turbulence near the leading edge of the blade results in an 
elevated level of k, strongly affecting the development of the turbulent boundary layer 
along the blade. Large flow turning and curvature effects, present in these leading edge 
flows, influence the development of the flow near the stagnation point. The experimental 
data on this effect, especially in turbomachinery stages, is scarce. The k-e turbulence 
model predicts the level of the turbulent kinetic energy. As a result, the boundary layer 
becomes fully turbulent, with the transition occurring very close to the leading edge. 

It is possible to separate the flow near the stagnation point into three regions. The 
first region is a freestream flow. In this region turbulent kinetic energy is balanced by the 
dissipation term in the k-equations. The second region is the thin boundary layer, 
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developing from the stagnation point. In this region the mean flow and the turbulence 
equations are strongly coupled and should be solved simultaneously. The third region is 
the buffer zone between previous two zones. Despite the fact that the mean flow can still 
be considered inviscid and the development of the turbulence does not affect the mean 
flow variables in this region, the velocity gradients severely affect the development of the 
turbulence. 

The modification of the k-e model to improve the accuracy have been suggested by 
many investigators. The first group of modifications suggested consists of change in the 
production term. The production of the turbulence kinetic energy can be expressed 
(incompressible flow): 


3U i dU: 2 dU - 

F - v, (i-r 1 - * T- 1 -) - t V>17 L - 2v > s ij • S H 


dx- dx j 


dXj 


[3-1] 


Launder (1974) suggested the use of the rotation rate to modify the production 
term: P = 2v l ^J(s ij JP ■ (r^ Y . The flow near (upstream) the stagnation point is nearly 


irrotational, while in the shear layer : -^(s^) 2 • (/? (J and: (s„) (s s ) are practically equal. 

This modification reduces the production of turbulence only near the stagnation point, 
while the boundary layer is not affected. Jin and Braza (1994) proposed the production 
term of the form: P = 2 v,R 9 ■ R tj . This modification can be used only in a limited region 
near the stagnation point However, the experimental data gives an increase of k up to 
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three times of the freestream value along the stagnation line. P = 2v l R IJ ■ R tj is zero along 

the stagnation line and correspondingly k will be about uniform. 

In the second approach, the dissipation rate equation is modified. The equation for 
the turbulence dissipation rate is one of the ’weakest’ parts of any turbulence model. The 
derivation of this equation is based on numerous assumptions. Strahle et al. (1987) 
suggested setting C e2 and Qi (Eq. 2-6, Chapter 2) near the stagnation point in the dissipation 
equation. This modification is based on the consideration of the analytical solution for k 
and e equations near the stagnation point. The potential solution was used to define 
velocities, i.e., it was assumed that the mean flow and turbulence equations are uncoupled 
in this region. According to Strahle et al. (1987), a consistent solution can be obtained 
only if C e2 and C £ i are equal. Due to the fact that Ce 2 is based on the decay of grid 
turbulence C £ i is chosen to be modified. 

Numerical simulations have been carried out to assess various modifications to the 
turbulence model indicated above. In the case of time marching scheme, stiffness of k-e 
equations and a highly stretched grid near the stagnation point lead to an additional 
difficulties near the leading edge. Viscous flow in the compressor cascade, described in 
Chapter 4, has been simulated using these modifications to the FLB k-£ model for the 
leading edge flow. The ratio of the normal to the shear stresses has been used as a switching 
function to switch from the modified Qi near the stagnation point to the original Qi in the 
shear layer. Results based on the original k-e model, without the modification, predict a very 
high level of the turbulent kinetic energy near the leading edge, while models with the modified 
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e-equation predict smaller production of k, closer to the observed values. The original model 
predicts very low turbulent dissipation rate near the stagnation point. Due to an excessive 
production of the turbulent kinetic energy and a low level of e, very large values of eddy 
viscosity are predicted. This leads to an excessive diffusion, which is not physical. As a result of 
this, the boundary layer becomes fully turbulent at the leading edge. All modifications improve 
the distribution of the turbulent kinetic energy and enable a reasonably good prediction of the 
transition. Predicted transition zone is located from about 33% of chord to 50% (this 
correlates well with the experimental data), with 5% variation between various 
modifications. The modification, due to Jin et al. (1994), leads to earliest transition 
inception while the modification of the e equation tends to predict the latest transition 
inception. Due to a lack of the experimental data on the budget of the turbulent kinetic 
energy along the stagnation line, it is difficult to assess various modifications suggested. 
Approaches, based on the modification of the e-equation and the modification of the 
production term, due to Launder and Spalding, have been chosen for the simulation of the 
unsteady transition described in following chapters and they are found to be crucial for the 
accurate prediction of the unsteady transitional flow. 

3.4.2 Turbulence flow field in the freestream 

The flow in turbine passage is another example of the ease with normal stress 
dominance. Rapid flow acceleration/deceleration outside boundary layers generate normal 
Reynolds stresses (in streamwise direction) that are higher than a shear stress. 


Comparisons between the predicted flow in Penn State turbine rotor (Chapter 5), based on 
k-e model, and the measured data shows that turbulence intensities may be overpredicted 
by more than ATu= 5% (Fig. 3-8). This problem has two main consequences. First is the 
effect on the boundary level transition. For the unsteady numerical simulation it also 
affects the decay of the upstream wake through the turbine passage. In contrast to the 
stagnation point flow, this problem is not contaminated by the near wall effect, thus results 
are independent of low-Re part of the model. Simulations based on the Fluent code have 
been used to get more insight into the turbulence flow at a midpassage of the turbine 
blade. Three different turbulence models are utilized for this analysis; standard k-e model, 
renormalized group k-e model (Orszag et al., 1993) and Full Reynolds stress model. 
Fluent tends to predict more radical rise in the turbulence level (ATu=2-3%) in 
comparison with Penn State code, therefore the comparison is done between different 
turbulence models using the same code (Fluent, Fig. 3-9 and Penn State, Fig. 3-8). 
Utilization of the RNG k-e model provides only a moderate improvement. This difference 
may be attributed to a modified coefficient of the RNG k-e model rather than to improved 
physics. Second-order turbulence model significantly improves the predicted distribution 
of the turbulent kinetic energy (Fig. 3-9). Error in the prediction based on the k-e model 
may be separated into two components; inadequate coefficient of the model and lack of 
the flow physics. A comparison of the streamwise and crossflow Reynolds stresses 
presented in Fig. 3-10 (note different scales for k-e and FRSM predictions), indicates that 
the predicted normal stress differs not only in amplitude, but has essentially different 
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distribution. Comparison of pseudo viscosity coefficient e at the location of the maximum 
turbulence intensity shows that e s ~ 2e„. Therefore the adjustment of the k-e model 
coefficient will not lead to the correct solution (e.g. RNG k-e model can only moderately 
improve the solution). Application of the second-order turbulence closure for 
turbomachinery flows is beyond the framework of this thesis. Hence, modifications similar 
to those analyzed for the stagnation flow have been considered for the k-e model 

modification. Only the modification of the prediction term based on P = 2 {r :j j 2 

is found to be suitable for the current task (Fig. 3-8). The modification with the prediction 
in the form: P = 2v r /?, y /? iy has been found to underpredict the upstream wake dissipation in 
the case of the rotor-stator interaction. 


tnf'l 


Fig. 3-1 Total velocity distribution along the stagnation streamline, inviscid simulation 
“H” type gird 



Fig. 3-2 Computational error 
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Fig. 3-5 Steady transition on a flat plate 
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Fig. 3-6 Steady transition on a flat plate, correlation from Abu-Ghannam and Shaw (1980) 
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Fig. 3-7 Effect of the artificial dissipation on the skin friction coefficient on a rotor blade 



Fig. 3-8 Turbulence intensity at midpitch inside the rotor blade passage 
data due Zaccaria and Lakshminarayana, 1997b 
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x/C, 

Fig. 3-9 Turbulence intensity at midpitch inside the rotor passage (FLUENT) 



Fig. 3-10 Predicted distribution of the normal Reynolds stresses along midpassgae 
streamline (FLUENT) 



Chapter 4 


NUMERICAL SIMULATION OF THE UNSTEADY TRANSITIONAL FLOW IN 

A COMPRESSOR CASCADE 


Simulations of the test cases, presented in the Chapter 3, establish the accuracy of 
the code for basic flows. Numerical analysis of the unsteady transitional flow in a 
compressor cascade, presented in this chapter, is carried out to assess the ability of the 
code to simulate realistic turbomachinery configurations. Even though attempts have been 
made to develop and use unsteady Navier-Stokes solvers for the prediction of rotor-stator 
interaction effects, none have been satisfactorily validated against accuracy, especially for 
its ability to capture the unsteady transitional viscous layers near blade and wall surfaces. 
This is the main objective of the research presented in this chapter. Characteristics of the 
different low-Reynolds k-e turbulence models have direct impact on the ability of the code 
to simulate the unsteady transitional boundary layer. Further research is required to assess 
the applicability of the k-e models for the unsteady flow modeling. Temporal accuracy of 
the code depends on the adequate choice of the number of inner cycles and the physical 
time iterations. Numerical simulation of the compressor flow is used to evaluate the 
influence of these factors on the accuracy of the prediction. 

4.1 Compressor cascade description 

Fan & Lakshminarayana (1996) used an unsteady inviscid two-dimensional code 
coupled with an unsteady boundary layer code to predict the unsteady flow caused by the 
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transport of a simulated upstream rotor wake through an annular compressor cascade and 
compared the prediction with the data due Schulz et al (1990). The simulated rotor wake 
was generated using a rotor consisting of rotating rods (Fig. 4-1). The annular cascade 
had twenty four untwisted blades. Characteristics of the cascade are given in Table 4-1. 
Fan and Lakshminarayana (1996) used 241(streamwise) X 61(blade-to-blade) grid for 
inviscid calculation. For the boundary layer solution, 120 streamwise stations and 121 grid 
point normal to the wall were used. 


Table 4-1 Compressor cascade characteristics 


Pitch/Chord 

0.78 

Stager angle 

29° 

Steady inflow angle, a 

44°, 49.2° 

Inlet Ma 

0.299 

Re 

4x 10 6 

Wake inflow angle 

15.55° 

Reduced frequency, Q 

6.12 

Wake width parameter, (0 

0.095 

Wake velocity defect, Ao 

28.3% 

Variation of the turbulence intensity, 
AT U 

8 % 

In the numerical 

simulation, tl 

ie inlet wake was prescribed as a 

Gaussian 


distribution, this was found to be a good approximation of the measured wake: 


V = Vq + A 0 exp( 2 ') 

2 co 

Similar distribution is used to prescribe an inlet distribution of turbulence characteristics. 
Amplitude and width are adjusted according to experimental parameters: 


(y-v ow t)\ 

T =T +AT exp( % ) 

u u0 u 2 co 1 


Where V ow - wake pitchwise velocity. 



72 


To assess the accuracy of the code, predictions based on the Navier-Stokes 
procedure are compared with the data and with predictions from the Euler/boundary layer 
approach. In some instances, the comparison is done only with the Euler/boundary layer 
prediction, due to a lack of the experimental data. Cases with 44° and 49° inlet flow 
angles have been chosen for numerical simulations. 

4.2 Sensitivity studies 

Flow simulations have been performed using three different grids: 179x61, 189x95 
(Fig. 4-2), and 201X193 to investigate the grid dependency. The distance between the first 
grid point and the wall varies from y+=1.6 for the coarse grid to y+=0.6 for the fine grid. 
Numerical simulation of the steady flow (C p and Q distributions) do not show any 

significant difference between predictions with 201x193 and 189x95 grids. Numerical 
simulation of unsteady flows impose additional requirements on the grid generation. Grid 
should be fine enough to allow a correct propagation and decay of the unsteady wake 
through the passage. The numerical analysis presented earlier (Chapter 3) is used to satisfy 
this requirement. Fourier decomposition of the inlet wake shows that it has five essential 
harmonics (based on blade passing frequency). Amplitude of the fifth harmonic is found to 
be only 1.3 % of an amplitude of the first harmonic. The grid with 193 grid points in the y- 
direction enables the wake to propagate through the cascade without non-physical decay, 
caused by numerical factors. In the case of 95 grid points in y-direction, only the fifth 
harmonic is affected by the artificial dissipation. This effect can be neglected, because the 
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fifth harmonic is dispersed rapidly by the physical dissipation. As a result of the grid 
dependency analysis, 189x95 grid has been chosen in all computations. 

The choice of the number of physical time steps per period and the number of inner 
iterations at each physical time step is an additional factor which affects the accuracy of 
the unsteady flow simulation. An increase in the number of physical time steps leads to a 
growth in the temporal accuracy and reduced phase errors. The main requirement is that 
the physical time step should be small enough to resolve the smallest time scales. Previous 
research indicated that about 500 physical time steps per period is required for the 
accurate solution of the wake-blade interaction effect. An increase in the number of 
physical time steps also affects the number of inner iterations, because of the smaller initial 
error. The change in the number of physical time steps modifies convergence characteristic 
of the scheme with a pseudo time stepping. For the case when At— **> (a steady state 
solution), the additional source term due to the presence of the physical time derivative 
vanishes. Additional damping at low wave numbers disappears and correspondingly, the 
convergence of the inner cycle in the freestream slows down. To analyze the effect of the 
number of physical and inner iterations on the accuracy of numerical results, a number of 
numerical tests have been carried out. 

Numerical simulations were performed with 500 physical time steps and 10, 20, 50 
inner iterations. About 1.5 order of magnitude drop in the maximum residual (mean flow 
equations) was achieved with 10 inner iteration at each physical time steps, this number 
increased to 2.5 for 50 inner iterations. A comparison of the predicted amplitude of the 
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first harmonic of the unsteady pressure (ACp) on the blade is shown in the Fig. 4-3. There 
is no significant difference between 20 and 50 inner iterations. Similar results were 
obtained for other harmonics and phase angles. On the other hand (Fig. 4-4), there is a 
significant change in the amplitude Cf variation on the suction surface, when the number of 
inner iterations is increased to 50. Additional numerical experiments showed that further 
increase in the number of inner iterations does not affect the accuracy of the unsteady 
C f distribution. Unsteady skin friction coefficient on the pressure surface is less sensitive 
to this factor because of the thin boundary layer. The convergence of the unsteady 
pressure depends on the convergence characteristics of the numerical scheme in the 
ffeestream. As it was indicated previously, inlet wake has five essential harmonics. This 
corresponds to the wave number range of 0.03 to 0.15. The ratio of the pseudo time step 
to the physical time step is high outside the boundary layer. As a result, an additional 
damping of the low wave number harmonics provides very rapid convergence of the 
unsteady ffeestream flow. During first 10 iterations the inner cycle convergence is equal to 
an analytical value based on the ffeestream Ax/At ratio. This fact also supports the 
previous conclusion. On the contrary, the correct prediction of the unsteady Cf requires an 
accurate simulation of the unsteady velocity in the boundary layer. In addition to 
employing smaller inner steps due to the fine grid, there is no ‘positive’ effect of the 
source term in this region (At/At - 0). Thus the prediction of the unsteady velocity 


requires more inner iterations. 
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Numerical simulations with 250, 500, 1 000 physical time steps were carried out to 
estimate the influence of the number of physical time steps. The number of inner iterations 
was 20 for all cases. According to the previous analysis, this number is sufficient to 
achieve a converged pressure field. The comparison, presented in Fig. 4-5, indicates that 
at least 500 steps are required for the accurate prediction. 

4.3 Unsteady pressure field 

The predicted pressure distribution and the unsteady pressure envelope on the 
blade are shown in Fig. 4-6. The time averaged blade pressure distribution from the 
Navier-Stokes code is more accurate than the prediction from the Euler code and is also in 
a good agreement with the experimental data. This is due to the presence of the separated 
region at about 90%-95% of the chord. The predicted time history of the unsteady 
pressure is compared with the experimental data in Fig. 4-7. There is a good agreement 
between predictions and the experimental data. Numerical simulation correctly predicts the 
maximum unsteadiness near the leading edge, where the wake hits the blade. This is 
caused by a change in the incidence angle and chopping of the wake by the blade. Beyond 
20% of the chord, the development of the recirculating flow pattern, induced by the 
passing wake, plays a dominant role in the development of the unsteady pressure. This 
leads to the region of an increased instantaneous pressure along the wake path. Fig. 4-7 
indicates that unsteady pressures are predicted well up to 20 percent of the chord, which is 
the most important and crucial part of the blade. Both the Euler and the Navier-Stokes 



76 


code predict smaller pressure variation in comparison with the experimental data, 
especially from x/Cx=60% to x/Cx=80% on the suction surface. This discrepancy is due to 
three-dimensional effects in the annular cascade. Flow visualization (Schulz et al, 1990) 
showed a strong comer separation. Interaction between the upstream wake and the comer 
separation leads to an amplification of the pressure oscillations. As expected, the Navier- 
Stokes solution predicts smaller amplitude in the unsteady pressure at the midchord in 
comparison with the Euler prediction due to the wake decay caused by the physical 
dissipation. There is a region of an increased unsteadiness in pressure near the trailing 
edge. This is due to the interaction between the passing wake and the separated region 
near the trailing edge. The Navier-Stokes code correctly predicts this feature (Fig. 4-5 and 
Fig. 4-7). While the flow is attached up to 97% of chord, velocity profiles from 85% of 
chord indicate a ‘near separation’ character of the flow. 

The development of the unsteady pressure field can be explained on the basis of 
two main mechanisms: wake cutting by the leading edge with the associated modification 
of the incidence angle and the development of the recirculating pattern due to the passing 
wake. Both phenomena are predominantly inviscid. As a result, both the Euler and the 
Navier-Stokes code predict nearly identical unsteady pressure field. 

4.4 Development of the unsteady transitional flow 

Prediction of the unsteady transitional flow is crucial in evaluating loses, efficiency 
and cooling requirements of turbomachinery. None of the Navier-Stokes procedures have 
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been validated for their ability to predict the unsteady transitional flow. Numerical 
simulations have been carried out using three low Re k-E models: CH, FLB and LB. The 
leading edge modification of k-£ models, described earlier, was found to be crucial for 
correct prediction of the transitional flow. Despite the modifications for the leading edge 
effect, the calculation with the CH model predicts fully turbulent flow all along the blade. 
Hence, the numerical simulation based on the CH model is not presented. 

As shown in Fig. 4-8a, in the boundary layer solution the transitional region 
extends from x/C x =0.1 to x/C x =0.30. An examination of the skin friction coefficient 
distribution in conjunction with the turbulence field shows that the Navier-Stokes code 
with the FLB model predicts transition from x/C x =0.12 to x/C x =0.4, while the 
computation based on the LB model predicts transition from x/C x =0.1 to x/C x =0.3. An 
understanding of the complex transitional process on the suction surface can be obtained 
from a study of the time history of the skin friction coefficient presented in Fig. 4-8. The 
trend predicted by both the Navier-Stokes code and the Euler/boundary layer code are in a 
very good agreement, unsteady fluctuations predicted by the Navier-Stokes code are 
slightly lower due to the decay of the wake, which is neglected in the Euler code. 

Wake induced transition is a very complex phenomenon driven by the interaction 
between the mean flow and the turbulence field. In this compressor cascade, the 
transitional region is located near the leading edge. Amplification and modification of the 
wake and turbulence due to the interaction with the leading edge affect the development 
of the transitional process. In Fig. 4-8b, path I corresponds to the upstream wake 
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propagation at the edge of the boundary layer, based on the maximum wake defect at that 
location. Path II (Fig. 4-8b) is the location of the maximum velocity fluctuations in the 
boundary layer. Beyond 20% of chord, an increase in phase lag between the convection 
velocity in the boundary layer and in the ffeestream is observed. An interesting 
phenomenon takes place in zone A (Fig. 4-8b ), along path II. At t/T=0, the wake is 
located at the leading edge (path I). Change in the incidence angle seriously affects the 
pressure distribution near the leading edge of the suction surface. As a result of this 
interaction, a zone with a reduced velocity is generated from 5% to 15% of the chord 
above the suction surface. This region modifies the development of the boundary layer. A 
zone of a low mean flow is located near the leading edge from about t/T=-0. 1 to t/T=0.15 
(Fig. 4-7 and Fig. 4-8) and disappears after passing of the wake. This phenomenon 
accounts for the difference in the location of the minimum Cf predicted by the Navier- 
Stokes and the boundary layer code shown in Fig. 4-8. Following the classification 
suggested by Halstead et.al. (1995), it is possible to identify various regions associated 
with the wake induced transitional flow. In Fig. 4-8b, A is the region of the wake induced 
transition. Disturbance due to the wake-boundary interaction leads to an earlier transition. 
The region B is a region with a transition between wakes and a zone with some features 
associated with the becalmed region. This region is located downstream of the small, fully 
laminar zone L. There is a smooth decrease in the shape factor . H in the zone B (line Z, 
Fig. 4-9), while in the region of the wake induced transition, a sharp drop of H (line Y) 
indicates an abrupt transition from the laminar to the turbulent flow. A comparison 
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between predictions based on the FLB model and the LB model shows that the LB model 
predicts an early transition (Fig. 4-9). Zone C is the wake induced turbulence strip. Two 
counterrotating vortices inside the boundary layer, generated by the interaction of the 
upstream wake with the cascade flow, have a major influence on the unsteady boundary 
layer downstream of the transitional region. A clockwise rotating vortex near the leading 
edge of the wake, leads to a smaller C f in the region C. The skin friction coefficient in the 
region D (trailing edge of the wake) has an increased level of Cf as a result of the 
counterclockwise vortex. Boundary layer simulation predicts no change in the skin friction 
coefficient in this region. The Navier-Stokes simulation provides more accurate simulation 
of the wake behavior in the outer region of the boundary layer and predicts the extended 
region of this counterclockwise vortex (about 20% of chord). Due to the presence of this 
vortex, the boundary layer profile has larger gradients near the wall resulting in increased 
shear stresses between wakes. 

The transitional flow and the unsteady boundary development are controlled by 
both the mean velocity defect and turbulence variation in the wake. These two factors 
have dissimilar influence in different regions. Unsteady interaction between the mean flow 
and the turbulence field, with a phase lag between the velocity, the pressure, and the 
turbulence quantities make the flow very complex. It is possible to estimate the 
importance of these two mechanisms from an analysis of the unsteady flow field simulated 
by the Navier-Stokes solver. In the wake induced transitional strip (zone A in Fig. 4-8 and 
Fig. 4-9), the influence of the velocity defect has the dominant influence. While the 
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turbulence intensity in the wake reaches 10%, the phase lag between k in the freestream 
and k in the boundary layer reduces the influence of the increased turbulent kinetic energy 
at this location. An opposite effect is felt in the region B, which is located between wake 
paths. The amplification of the unsteady turbulence near the leading edge leads to an 
increased level of k from 10% to 30% of the chord. Additional diffusion of the turbulence 
from the freestream results in a smoother transition, as seen along line Z in Fig. 4-9b. 

One of the characteristic features of the unsteady boundary layer is the phase lag 
between the velocity at the edge and inside the boundary layer. The distance between path 
I and path II (Fig. 4-8b) is widening with the development of the boundary layer 
downstream of the leading edge. This is an indication of an increased phase lag in the 
velocity field. The amplitude and the phase angle of velocity fluctuations are shown in Fig. 
4-10. This phase lag increases from about 30° at Xc/C=0.4 to 100° at Xc/C=0.76 (Fig. 4- 
10b). The predicted phase angle and the amplitude of the first harmonic of the total 
velocity correlate well with predictions from the boundary layer code. The Navier-Stokes 
solution predicts sudden increase of the phase lag in the laminar sublayer. This can be 
attributed to an inadequate grid resolution in this region. It should be remarked here that 
the boundary layer code employs 121 grid points inside the boundary; nearly four times as 
many as that used in the Navier-Stokes solver. Considering this, the agreement is good. 
This is one of the most important and critical steps in the validation of the Navier-Stokes 
code. With proper control of the artificial dissipation, grid, and time step, the Navier- 
Stokes code can be used to predict the unsteady transitional boundary layer accurately. 
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The predicted momentum thickness, shown in Fig. 4-11, reveals an excellent 
agreement with the experimental data and the Euler/boundary layer prediction. Similar to 
the boundary layer solution, the Navier-Stokes solver predicts a higher level of the time- 
averaged momentum thickness in a comparison with the steady state solution. This is an 
indication of an increased loss due to the unsteady interaction. According to the flow 
visualization, the flow separation occurs near 90% chord. The main drawback of the 
boundary layer approach is its inability to simulate the separated flow. The Navier-Stokes 
code correctly predicts separation zone, existing from 87% of chord. The predicted 
separation has an unsteady character, flow conditions vary from the fully attached to the 
separated flow. A sharp increase in the momentum thickness beyond 85% of chord is due 
to an earlier separation caused by the passing wake. 

4.5 Stator wake 

The development of the stator wake is influenced by its interaction with the 
upstream rotor wake. As described previously, the passing wake generates two 
counterrotating vortices which have considerable influence on the development of both the 
pressure and suction side boundary layers. The amplitude and extent of the secondary 
vortices inside the pressure surface boundary layer are smaller due to extremely thin 
boundary layers. Vortices inside the suction surface boundary layer play a dominant role in 
the development of the unsteady wake. This is evident from the unsteady velocity and 
turbulent kinetic energy distribution in the static wake shown in Fig. 4-12, Fig. 4-14, and 
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Fig. 4-13 respectively. The unsteadiness is higher away from the wake center. This 
distribution is very similar to that inside the boundary layer shown in Fig. 4- 10a. The 
fluctuations are zero at the wall and this is reflected in a very low unsteadiness inside the 
wake center. Maximum unsteadiness in kinetic energy occurs on the suction side of the 
wake and the reasons for this are explained below. 

The convective speed of the upstream wake is different for the pressure and 
suction surfaces. At the trailing edge, the phase angle difference between the passing wake 
near the suction surface and the pressure surface approaches 100°. This is an additional 
source of the unsteadiness in the wake. Clockwise vortex generated inside the boundary 
layer due to the wake-boundary layer interaction is shed into the stator wake, thus 
amplifying unsteadiness due to the wake passing. Clockwise rotating vortex, located 
above the stator wake, produces a system of counterrotating vortices in the stator wake. 
The presence of the unsteady vortices in the boundary layer creates an additional unsteady 
amplification of the wake turbulence. Fig. 4-13 indicates that this amplification is much 
higher on the suction side of the profile wake in comparison with the pressure side. The 
unsteady total velocity and turbulence kinetic energy distribution in the wake (Fig. 4-12 
and Fig. 4-13) indicate that the disturbance due to the passing wake is mainly concentrated 
in the first harmonic. The effect of the counterrotating vortices and its shedding is mainly 
confined to second and third harmonics. Interaction between vortices also leads to an 
amplification of the spatial oscillation of the wake. The deviation in the trajectory of the 


center line of the wake is about 3% of the chord at x/Cx= 1.29. 
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The interaction between the passing wake and the profile wake results in complex 
vortex structure, increased unsteadiness in the profile wake and high level of turbulent 
kinetic energy. These effects result in increased dissipation and faster decay of the stator 
wake. Thus the presence of the unsteadiness due to the rotor-stator interaction tends to 
decay profile wake faster. The time history of the instantaneous velocity, shown in Fig. 4- 
14, clearly reveals the location and the extent of the maximum interaction region. This 
occurs when the passing wake is directly inside the profile wake. The effect of additional 
shed vortices, causing additional oscillations, is also evident from this plot. 
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Fig. 4-1 Scheme of the experiment 



Fig. 4-2 Computational grid 
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Fig. 4-3 Influence of number of inner iterations on the amplitude of unsteady pressure 



Fig. 4-4 Influence of number of inner iterations on the amplitude of skin friction coefficient 












a) Navier-Stokes prediction 



0.0 0.2 0.4 x/Cx °' 6 08 1 *° 



£ 


c) Experiment 



0.0 0.2 0.4 x/Cx °' 6 °‘ 8 


Cp-Cp aver 

0.10 

— 0.08 
— 0.06 

R 0 04 


1.0 


Fig. 4-7 Unsteady pressure distribution, suction surface 
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Fig. 4-8 Skin friction coefficient 
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b) Navier-Stokes prediction FLB model 
Fig. 4-9 Shape factor 
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Fig. 4-1 1 Momentum thickness, suction surface 
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Fig. 4-12 Amplitude of the total velocity in the stator wake 






Chapter 5 


UNSTEADY FLOW IN TURBINE 

Rotor-stator interaction in a turbine stage has more complex structure in 
comparison with the unsteady flow in the compressor cascade analyzed in the previous 
chapter. In a compressor, the upstream wake propagates as a set of straight segments. The 
shape of each wake segment undergoes only minor distortion when the wake is 
transported through the compressor stage. Significantly higher turning and acceleration of 
the flow in a turbine stage results in a strong distorution, twisting and stretching of the 
wake segment. There is an increased phase lag between the flow at the midpitch and at the 
edge of the boundary layer. The difference in the relative wake-flow alignment also 
contributes to the modified rotor-stator interaction in a turbine stage. The ultimate 
solution of the rotor-stator interaction has to come from a systematic, scientific, and 
building block approach. The measurement in an actual engine is not only complicated, but 
will rarely provide an insight into numerous sources of unsteadiness. Likewise, a 
computational code with artificial dissipation and numerical error may mask some of the 
important physics. The major objective of this research is to understand the nature and 
magnitude of unsteadiness in turbines and develop insight for incorporation of these 
effects in the design and analysis. This objective is accomplished through a coupled 
experimental and computational study. Author has carried out the computational part and 
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participated in the interpretation of the experimental data. The experimental data was 
acquired by D.Ristic (Lashminarayana et al., 1998). 

The coupled experimental and computation study carried out at Penn State, which 
is more comprehensive than hither to attempted, should provide a detailed understanding 
of the flow physics. Almost all earlier experimentalists have concentrated their effort on 
blade surface measurements (e.g., pressure, skin friction, and transition). The present 
approach is to derive information on both the velocity field and the blade surface 
properties due to rotor-stator interaction. Likewise, the computational approach involves 
a development of a two-dimensional unsteady Navier-Stokes code, incorporating the 
transition and turbulence models to predict the nozzle wake, unsteady rotor blade 
boundary layers, unsteady free stream flow field and blade pressures accurately, and 
carrying out a simulation study to understand the sources of unsteady viscous flow 
through turbomachinery stages, including its effects on the transition and aero- 
thermodynamic losses. 

5.1 Experimental Program 

The Axial Flow Turbine Research Facility (AFTRF) of The Pennsylvania State 
University is an open circuit facility 0.9166 m (3 feet) in diameter and a hub-to-tip radius 
ratio of 0.73, with an advanced axial turbine blading configuration (Fig. 5-1). The facility 
consists of a large bellmouth inlet, followed by a test section with a nozzle vane guide row 
and a rotor. There are 23 nozzle guide vanes and 29 rotor blades followed by outlet guide 
vanes. A window for LDV measurements covering the entire flow field from upstream of 
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the nozzle to downstream of the rotor passage is also incorporated. Detailed design of the 
facility, performance, and geometric features are described in Lakshminarayana et al. 
(1992). Some important performance and geometric parameters are as follows: ratio = 
0.7269; Nozzle, axial chord (midspan) = 11.23 cm, turning angle = 70 degrees; Rotor : 
hub/tip axial blade chord (midspan) = 9.294 cm, turning angle = 95.42 degrees at tip, and 
125.69 degrees at root, tip clearance = 0.97 mm; Reynolds number of nozzle flow (based 

on exit flow) = 10 6 , mass flow rate = 11.05 kg/s, loading coefficient (2A P a I pU^) - 3.88, 

rotational speed = 1330 rpm. The vane-blade spacing is 22.6 % of the nozzle axial chord 
at the midspan. The design velocity triangles at the inlet and exit of the rotor are shown in 
Fig. 5-1. 

Comprehensive data was acquired within and downstream of the nozzle during the 
earlier phase of the research program. A complete five-hole probe survey was carried out 
at 2.5 and 9% nozzle chord downstream of the nozzle trailing edge (Zaccaria and 
Lakshminarayana, 1995). Furthermore, LDV measurements carried out at midspan of the 
rotor (upstream to downstream) establishe the details of the distortion and the nozzle 
wake profile upstream of the rotor (Zaccaria and Lakshminarayana, 1997b). Hence all the 
rotor inflow properties (time dependent in rotor frame) are known and shown in Fig. 5-2. 
The LDV measurements at the rotor midspan are carried out at the design condition 
(Table 5-1). The blade dynamic pressure data at midspan and hub are acquired at design 
and three off-design conditions, and these are shown in Table 5-1. This was achieved by 
varying the turbine speed and the mass flow. The off-design data listed in Table 5-1 is 



96 


based on the one-dimensional consideration. The upstream flow field was not acquired at 
the off-design conditions. The inlet flow properties (x/C x = -0.224) are shown in Fig. 5-2 
at mid span (H = 0.5) of the rotor and near the hub (H = 0.035). These are also the inlet 
conditions used for the computation. It is clear that substantial velocity, flow angle, 
stagnation, and static pressure gradients exist in the nozzle wake, both at the hub and the 
midspan regions. Most of the earlier research was carried out at larger nozzle-rotor 
spacing, where the pressure wakes do not exist as they decay rapidly downstream. The 
present measurement is carried out at practical spacing; and hence, all the features of 
modem turbines are present, including both the pressure and velocity gusts due to the 
viscous wakes. The flow field near the hub includes secondary flow regions as evidenced 
by the defect in velocity, static pressure, and stagnation pressure near the outer edge 
(suction side) of the wake. 

Table 5-1 Operating Conditions 


Operating Conditions 

mm 

B 

c 

D 

Mass flow (Kg/s) 

10.42 

10.42 

5.82 

5.82 

RPM 

1330 

1235 

1100 

800 

V x (ms) 

30.43 

30.43 

16.8 

16.8 

U m (ms) 

55.32 

51.35 

45.75 

33.27 

Rotor inlet flow angle , Pi , (deg) 

41.44 

45.76 

-4.75 

33.4 

Incidence (deg) 

-1.9 

2.06 

-44.9 

-7.6 

£2 = GJ C /V 

X x 

4.89 

4.54 

7.32 

5.82 

Re 



ilium 

iimiijfl 


The reduced frequency based on the upstream nozzle wake frequency, rotor axial 


chord and axial velocity upstream of nozzle varies from 4.54 to 7.32. Since there are 
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significant differences between the condition A (design) and C (off design), only the data 
at these two conditions are presented and interpreted. 


5.1.1 Laser Doppler Velocimeter Data 

Detailed steady and unsteady flow data were acquired at midspan of the rotor 
blade using a two-dimensional laser doppler velocimeter. Since the flow is two- 
dimensional at this spanwise location, only the axial and tangential velocity components 
were measured. The measurements were carried out at 37 axial locations upstream of the 
rotor (x/C x =-0.088) to one chord downstream of the rotor (Zaccaria and 
Lakshminarayana, 1997b). To account for the non-uniformity of the rotor absolute inlet 
flow field, measurements were made at six tangential locations in the absolute frame 
equally spaced over one nozzle pitch (Fig. 5-1). These six tangential positions represent 
six different relative positions between the nozzle and the rotor (labeled position 1 through 
6) or if viewed from the nozzle frame of reference, six different time-resolved positions of 
the rotor in relation to the nozzle. 

At each survey point, approximately 120,000 velocity measurements were 
acquired. Since all the velocity components were spatially phase-lock averaged, which 
result in a representative rotor passage with 50 measurement windows, there were 2400 
velocity measurements on average in each measurement window. After all instantaneous 
velocity measurements were acquired for each particular survey point, the velocity was 
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ensemble-averaged at each measurement window. The unresolved velocity for each 
measurement window can be calculated as 


v' = V i -V [Eq. 5-1] 

where V) is the instantaneous velocity measured at a particular rotor measurement 
window. The variance is given by: 


(v'Y= L 


n _ 

X (V. -V) 2 
1 = 1 1 


(fl-1) 


[Eq. 5-2] 


The level of unresolved unsteadiness in each measurement window is determined 
by the variance. For the LDV measurements in the rotor, the instantaneous velocity, V in , 
is decomposed as follows 

V in = V+V+V' [Eq. 5-3] 

where V is the time-averaged velocity, V is the periodic velocity and V’ is the unresolved 
velocity component as calculated in Eq.5-1 . 

The cycle-averaged values were obtained by averaging the ensemble-averaged (and 
phase-lock averaged) properties in each rotor measurement window for one nozzle/rotor 
location over the six nozzle/rotor locations (see Fig. 5-1). 

It should be emphasized here that the experimental procedure and data processing 
are designed to obtain spatially and “temporal” measurements (i.e., rotor shaft positions; 
not real time) of the wake-rotor interaction generated unsteadiness in the rotor. The laser 
is located at a fixed position relative to the nozzle wake for each nozzle/rotor location. 
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Hence, the averaging is based on identical rotor-nozzle blade positions. These 
measurements are similar to those reported by earlier investigators (e.g., Hathaway, 
1986). 

The LDV measurements are subject to numerous errors, most of which can be 
quantified. A complete error analysis for these measurements is given in Zaccaria and 
Lakshminarayana (1997a). Based on this error analysis, the uncertainty for a 95% 
confidence level is as follows: outside the rotor wake 0.4% and 2.8%, respectively for the 
ensemble-averaged velocity and the unresolved component of velocity; inside the rotor 
wake, 4.0% and 14.8% for the ensemble-averaged and unresolved velocity, respectively. 
The LDV data is presented and interpreted in Zaccaria and Lakshminarayana (1997b). 
This data is used in this thesis to validate the code and to derive new insight on the nozzle- 
rotor viscous interaction in view of the additional blade surface data acquired recently and 
the comprehensive flow simulation carried out with an unsteady Navier-Stokes code. 

5.1.2 Dynamic Pressure Measurements 

To provide a complete knowledge of the nozzle-rotor interaction phenomena the 
instanteneous unsteady dynamic pressure has been acquired (Lakshminarayana et al., 
1998). The sensitivity and the high-frequency response of silicon based semiconductor 
strain gauge pressure transducers coupled with their small size make them very desirable 
for use in obtaining unsteady dynamic pressure fluctuations. The implementation of the 
dynamic pressure transducers in the turbine rig was driven by space limitations. The 
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pressure transducers are inserted into chambers, which in turn are connected to the turbine 
airfoil surfaces through 0.5 - 0.8 mm diameter holes. The design objective is to achieve a 
frequency response of 40 kHz. The miniature sensors used are the Kulite model XCS-093 
with a pressure range of maximum 5 psia. They are capable of measuring pressure 
fluctuations to an accuracy of 0.01 psia. Sixteen Kulite transducers are located along the 
chord at the midspan of the rotor blade. Seven transducers are on the pressure side and the 
remaining nine are located on the suction surface of the next blade in the same passage. A 
schematic showing the location of dynamic pressure transducers on the blade at midspan 
(16) and on the hub (endwall) surface (5) is shown in Fig. 5-1. The low-level signals from 
the dynamic pressure transducers are amplified in the rotating frame by using miniature 
amplifiers. The amplifiers rotate in the rotor frame and provide a high-level signal output 
before the signal reaches the slip-ring unit. These amplifiers are located inside the rotating 
instrumentation drum. The transducers were calibrated by inserting the entire blade in a 
pressure chamber as well as using the steady state static pressures on the blade before the 
experiment. 

The dynamic pressure data from the rotor is transmitted through the rotating drum 
to a slip-ring unit. The slip-ring unit is of the brush type and has 150 channels. The slip- 
ring unit is housed in a cowling in front of the facility. Each ring carries four brushes made 
of silver graphite. The rings are made up of coin silver, which .withstands up to current 
levels of 5 Amps. The brushes are individually removable and replaceable. The contact 


resistance is about 5 milliohms maximum. 
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A completely automated data processing system is built around a computer with a 
clock rate of 7MHz. The data is triggered by an encoder on the turbine shaft. 
Approximately 167 (1 million samples at 6000 samples/rev.) revolutions of data are 
acquired. This was ensemble-averaged to derive revolution periodic, which is then 
decomposed into blade periodic (average-passage) and blade aperiodic pressures. The 
procedure used for data processing is similar to that described in the last section. The 
latter represents the deviation of the individual passage from the average passage. A 
typical set of processed data for x/C x =0.285 is shown in Fig. 5-3. The blade aperiodic 
quantity is small for most passages and appreciable only in three or four passages. Only 
the blade periodic data are presented and interpreted here. 

5.2 Numerical Procedure 

To simulate rotor unsteady flow, time dependent boundary condition were 
imposed at rotor inlet. Velocity and the pressure distributions were based on the 
experimental data at 24.59% of the rotor chord upstream of the leading edge. Data was 
obtained using a five-hole probe in a stationary frame of reference. To impose an accurate 
unsteady boundary condition, the measured pressure and velocity fields (moving in rotor 
frame of reference) were superposed with the rotor potential field based on a steady 
prediction (steady in the rotor frame of reference). 

In order to apply periodic boundary conditions, the flow in five rotor passages was 
calculated. This provides a vane-to-blade ratio equal to 4:5, close to the actual ratio of 


23:29. 
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5.3 Sensitivity Studies 

To analyze the grid sensitivity, numerical simulations have been carried out using 
three grids: 101X61, 181X91, and 211X181. The distance between the first grid point and 
the wall varied from y+=0.9 for the coarse grid to y+=0.3 for the fine grid. Steady flow 
predictions based on 181X91 and 21 1X181 girds are identical. A numerical simulation of 
the unsteady flow imposes an additional requirement on the grid density outside the 
boundary layer. An earlier analysis of the influence of the artificial dissipation (Chapter 3) 
has been used to estimate the required mesh distribution in the middle of the passage. The 
axial gap in the turbine stage is equal to 27.6% of the rotor chord. Hence, the nozzle 
wake defects are substantial. At least ten Fourier harmonics are needed for an accurate 
representation of the inlet velocity field. The amplitude of the 10th harmonic is 2% of the 
wake defect. According to the previous analysis, the finest grid (211X181) provides 
enough spatial resolution to model the whole spectrum of the inlet wake. A numerical 
simulation based on 181X91 grid may lead to excessive decay of the 8th and higher 
harmonics due to the artificial dissipation. A comparison of the numerical prediction of the 
unsteady flow based on these two grids indicates that utilization of the coarser grid 
(181X91) results in a 3% smaller local wake defect due to the artificial dissipation. No 
other significant differences have been observed. To minimize CPU time, most of the 
simulation studies were carried out using 181 X 91 gird. 

The accuracy of the numerical simulation depends on the appropriate choice of the 
number of physical time steps as well as the number of inner pseudo-time steps. Current 
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simulations were carried out with 1000 physical time steps per period. Even though it was 
found that 500 time steps were enough for the accurate temporal resolution, the solution 
had features of numerical instability. This problem was eliminated when the number of 
physical time steps was increased to 1000. An investigation of the dependency of the 
solution on the number of inner iterations indicated that an accurate prediction of the 
unsteady pressure field required about 7 inner iterations. An accurate resolution of the 
velocity and turbulence field required 15 inner iterations. 

5.4 Time-Averaged Flow Field 

A comparison between the predicted (steady) and measured surface pressure 
distribution is shown in Fig. 5-4 and Fig. 5-5. The comparison is done for four cases 
(summarized in Table 5-1). Very good agreement between the predicted and the 
experimental data is achieved for the design point (Fig. 5-4) as well as for the flows with 
moderate incidence angle (Fig. 5-5). For the case C with the inlet flow angle equal to 
-4.75°, the numerical simulation predicts flat pressure distribution along 60% of the 
chord on the pressure surface. The experimental data indicate that the pressure on the 
pressure surface decreases up to 30% of the chord and increases from 30 to 60% chord. 
This discrepancy is attributed to the known deficiency of the k-e turbulence model in 
accurately predicting the extent of flow separation. The three-dimensional nature of the 
separation bubble may also contribute to this discrepancy. At ^=-4.75°, the predicted 
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separation bubble extends from x/C x =0.1 to x/C x =0.6. Its thickness is about 25% of the 
blade spacing at x/C x =0.3. 

A comparison of the blade surface pressure distribution at different operating 
conditions indicates that the flow is not sensitive to the moderate variation in the incidence 
angle and the prediction is excellent at these conditions. A lower inlet flow angle moves 
the stagnation point towards the suction surface. For the largest incidence angle, the 
stagnation point is located at x/C x = 2%. At P =-4.75° the flow near the leading edge 
undergoes significant modification. A strong separation bubble develops on the pressure 
side of the blade due to high turning of the fluid particles near the leading edge. The 
resulting pressure field leads to a compressor-like behavior from the leading edge to x/C x 
= 18% of the chord. 

The predicted (time-averaged) and measured (cycle-averaged) relative total 
velocity distribution at several axial locations at design condition (A) are shown in Fig. 5- 
6. The prediction is excellent at all locations except at x/C x = 1 10%. One of the severest 
tests on the code is the ability to capture a large gradient in velocity. Near the leading 
edge, the excellent agreement at x/C x =-5% provides confidence in the code. The wake 
depth is captured at x/C x = 110%, but the predicted wake width is smaller than the 
measured one. This may be attributed to three-dimensional effects or the cycle averaging 
procedure. The predicted passage average angle at the exit is 67.4°. This compares well 
with the measured and design values of 64 . 3 ° and 67 ° , respectively. These results indicate 
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that the time-averaged flow field is predicted accurately especially in the passage and 
upstream of the passage. 

5.5 Unsteady Pressure Field on the Blade Surface (Design) 

Before commenting on the nature and magnitude of unsteady pressure, it should be 
remarked here that the blading for the Penn State HP turbine was designed to limit over- 
acceleration on the suction surface just off of the leading edge (Lakshminarayana et al., 
1992). On any turbine blade, the flow is subject to high acceleration around the leading 
edge region on both the suction and the pressure sides. In both cases, there may be a 
subsequent undesirable diffusion. However, a local diffusion, however small, is of more 
concern on the suction side where it may constitute a significant disturbance to the 
(laminar) boundary layer. In the present blading design, it was possible to trade one off 
against the other by precise tailoring of the leading edge geometry and adjustment to the 
leading edge inclination (i.e., incidence). Generally, only a very slight tendency to 
overspeed was tolerated on the suction surface as it was felt that a very minor increase in 
the pressure surface diffusion would not significantly affect the blade performance. The 
suction peak is located near the mid-chord location. 

The measured and predicted steady pressure and unsteady envelope at the design 
condition as well as steady pressures are plotted in Fig. 5-4. The predicted amplitude of 
unsteady pressures agrees quite well with the measured values. The amplitudes are small 
near the leading edge on both surfaces. This is due to the fact that the blading is not very 
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sensitive to incidence changes due to inherent design. The incidence angle changes (Table 
5-1, Fig. 5-5) from 2.06° to - 7.6° had no significant input on the steady pressures. The 
maximum amplitude occurs near 28% chord on the suction surfaces and decreases 
gradually to insignificant values as the trailing edge is approached (Fig. 5-4). There is a 
significant difference between this data and those due to Dring et al. (1982), who observed 
maximum values at the leading edge. This may be due to large flow acceleration near the 
leading edge and closer rotor-stator spacing. The amplitudes (Fig. 5-4) are small on the 
pressure surface and once again decrease gradually to small values as the trailing edge is 
approached. 

It should be remarked that even though the agreement between the measured and 
predicted steady pressures on the suction surface is not exact from 30% to 50% chord, the 
amplitudes of unsteady pressure are in good agreement. 

The periodic variation of blade pressures on the suction surface at selected rotor 
chord locations are shown and compared with the predictions in Fig. 5-7. As mentioned 
earlier, the amplitudes are highest near the x/C x =0.28 and decrease to insignificant values 
near the trailing edge. The prediction is very accurate at x/Cx=0.071 and 0.28. 
Discrepancies are observed aft of the midchord. The measurement shows larger decay in 
the amplitude than the prediction. 

The space-time distribution of unsteady pressures (A C p ) on both surfaces along 

all measurement locations are shown in Fig. 5-8. Such plots can provide a valuable 
overview of the flow physics associated with the wake-rotor interactions. On the suction 
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surface, the measurement shows a maximum value at about 28% of the chord. The phase 
angle of the peak values (marked HS) on the suction side changes rapidly as the flow 
progresses downstream due to large changes in the convection velocity of the nozzle 
wake. The first low value on the suction side (LSI) also moves at about the same phase 
as the peak (HS). But the second minimum (LS2) is nearly at constant phase. On the 
other hand, the peak values on the pressure side occur near the leading edge. One of the 
unusual features is the presence of two peaks (HP1, HP2) on the pressure side. As 
indicated earlier and in Fig. 5-2, the inflow has both the pressure gust and the velocity 
gust. As the two approach the leading edge, the phase angle between the two distortions 
change and the pressure gust seems to have a more profound effect on the pressure 
surface than on the suction side. Another interesting feature is that, except for the leading 
edge, the peaks move at nearly constant phase as the flow progresses downstream along 
the pressure surface. Some of these features are similar to those observed by Dring et al. 
(1982), but major differences are in the location of peaks on the suction side. This is due 
to the differences in the blading design. 

The first three harmonics of the unsteady pressure on the suction surface are 
plotted and compared with the predictions in Fig. 5-9. The experimental data indicate that 
the first harmonic is the dominant component and this component is predicted quite well 
by the code. The predictions reveal two maxima; one near the leading edge and the other 
near 28% chord downstream. Since there is no data available at the leading edge, the 
presence of the first peak cannot be confirmed. But the minima observed near the 5% of 
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the chord location is predicted accurately. The second harmonic shows a peak close to the 
50% of the chord, but the predicted peak is further downstream. The occurrence of the 
peak in second harmonic can be attributed to the transfer of energy from the first harmonic 
to the second harmonic in the freestream unsteadiness. This is clear from the fact that the 
first harmonic decreases rapidly near the midchord and the peak in the second harmonic 
appears at this location. This is due to the dominant effect of two counter-rotating 
vortices (see discussion below) that are developed on either side of the wake and the 
weaker nozzle wake effect. This effect is overpredicted by the code. During the numerical 
simulations, it was observed that the amplitude of the second harmonic is sensitive to the 
local wake velocity defect. The measured values of the third harmonic are negligibly 
small. 

5.6 Unsteady Pressures at Off-Design Conditions 

The space-time history of unsteady pressures at the off-design condition (Case C, 
Table 5-1) is shown in Fig. 5-10. The relative flow at this operating condition is nearly 
axial, and the incidence is - 44.9° (Fig. 5-1). The peak unsteady pressures on both 
surfaces occur near the leading edge. One of the salient features of the distribution on the 
suction side is the change in the phase angle of the unsteady pressure as flow propagates 
downstream. The blade up to 20% of the chord behaves like a compressor. The relative 
velocities on the suction side of the turbine (pressure surface) are small up to about 20- 
30% of the chord (Fig. 5-5, Case C). Hence, the disturbance is nearly at constant phase 



109 


up to about 30% of the chord on the suction side. Significant changes occur beyond this 
location and the unsteady pressure decays very rapidly. This rapid decay is caused by 
major viscous effects within the passage due to the severe off-design conditions. The 
unsteady pressures decays more rapidly on the pressure surface. This is caused by a 
possible flow separation and a bubble near the leading edge of the pressure surface at this 
operating condition. The unsteady pressure on the pressure side (such as the suction 
surface) increases significantly. 

The measured values are compared with predictions at selected chordwise 
locations in Fig. 5-11. Here again, the agreement is very good near the leading edge, 
especially on the suction surface. As indicated earlier, the predicted flow separation on 
the pressure side may not be accurate, and this may account for the discrepancy at 
x/C x =0.3 on the pressure surface. 

The chordwise distributions of the unsteady pressure amplitude at all operating 
conditions (Table 5-1, Cases A, B, C, and D) are shown in Fig. 5-12. There are two 
major effects involved here; the time-averaged incidence changes and the reduced 
frequency changes. Cases A and B are carried out at the same speed, with Case B being 
at a higher mean incidence. The data on the suction surface shows that the amplitudes are 
slightly higher for Case B at most locations. The dominant effect is downstream of the 
leading edge (x/Cx=0.28 on suction surface, 0. 14 on the pressure, surface). Cases C and D 
have the same mass flow and Reynolds number, with Case C at more severe off-design 
condition. The data for Case D follows trends very similar to those of Cases A and B. As 
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indicated earlier, due to large negative incidence for Case C, the peak unsteady pressures 
are similar on both the suction and the pressure surfaces and have the fastest decay on the 
pressure side. The amplitudes are highest for this case on the pressure side. One of the 
interesting features of all of this data is that the peak pressure for all operating conditions 
(Cases A, B, C, and D) occur downstream of the leading edge on both surfaces. This has 
a significant impact on the leading edge cooling. The present design seems to have 
reduced unsteady wake-interaction effects near the leading edge through aerodynamic 
design. 

5.7 Unsteady Flow Field Due to Rotor-Stator Interaction 

A brief description of the data acquired from the LDV inside the passage at the 
midspan of the turbine rotor is given above. This data is compared with the predictions in 
this section. 

To account for the non-uniformity of the rotor absolute flow field, measurements 
were made at six tangential locations, equally spaced over a nozzle pitch (Fig. 5-1), in the 
absolute frame of reference. These six tangential locations represent six different relative 
positions between the nozzle and the rotor. The derived flow field corresponds to the 
fixed “rotor shaft positions,” not real time. More details on the experimental procedure 
and experimental measurements can be found in Zaccaria and Lakshminarayana (1997b). 
To enable a comparison of the measured “rotor shaft phase-locked” blade-to-blade flow 
field with the numerical prediction, the data from the computation has been processed in 
the same way as the experimental technique. The derived blade-to-blade distributions of 
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the wake defect and the unresolved unsteadiness (eq. 5-2) in the absolute frame of 
reference are shown in Fig. 5-13 and Fig. 5-14, respectively. It should be remarked here 
that the unresolved unsteadiness consists of both the random fluctuations due to 
turbulence as well as the blade periodic fluctuations, which is found to be small in this 
single-stage turbine. The velocity field is shown for the position N1 (Fig. 5-13, bottom 
station); while the unresolved unsteadiness distribution is plotted for the position N4 (Fig. 
5-14, top station). Due to the blade shadow and lack of seeding, no experimental data is 
available in the unshaded (white) region. The velocity represents the periodic unsteadiness 
due to the nozzle wake V (ensemble average minus the time average, eq. 5-3). A 
comparison between the data and the prediction shows good agreement with the location 
and defect in the nozzle wake. The numerical prediction moderately under-predicts the 
level of the wake defect between x/C x =0 and x/C x = 0.2 (Fig. 5-13, Fig. 5-15). The 
measured unresolved unsteadiness and the predicted turbulence intensity, (Fig. 5-14) 
indicate that the peak intensities are predicted reasonably well, but the wake width based 
on the unsteadiness shows that the computation has a larger diffusion (into the ffeestream) 
compared with the experimental data. 

Even though the “rotor shaft phase-locked” flow field can be used to analyze the 
development of the unsteady flow in a rotor passage due to the wake-rotor interaction, it 
is more appropriate to consider the instantaneous flow field. “Rotor shaft phase-locked” 
representation is used only to compare the measured velocity field with the numerical 
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prediction. All subsequent discussion is based on the time accurate (instantaneous) flow 
field. 

5.8 Discussion of Unsteady Flow Physics 

As indicated earlier, the vane-blade ratio used in the computation is 4:5; hence the 
unsteady flow field is computed in five rotor passages. Hence, all the simulation data 
shown in this and subsequent sections shows predictions in five rotor passages. 

The unsteady velocity, pressure, turbulence field, and the unsteady relative total 
pressure coefficient shown in Fig. 5-16, Fig. 5-17, Fig. 5-18, and Fig. 5-19 respectively, 
reveal the development of the rotor unsteady field caused by the nozzle wake-rotor 
interaction. Based on the unsteady velocity, the nozzle wake appears as a negative jet 
moving towards the suction side (Fig. 5-16). The nature of such impingement depends on 
the inlet velocity triangle and the wake defect. The location and propagation of the nozzle 
wake can be clearly identified by examining all the instantaneous properties: unsteady 
velocity (V -V), static and stagnation pressure, and turbulence field. The shape of the 
turbulence wake is very close to that of the entropy wake, which is normally used to 
identify the wake position. For the sake of brevity, the unsteady entropy distribution is not 
shown here. 

The wake segment is initially straight (Fig. 5-16). Due to the potential effect and 
the blockage caused by the leading edge, it is stretched, distorted, and diffused as it 
approaches the rotor leading edge. Further downstream, the wake is chopped and 
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transported at the local convection speed as a separate segment. Due to a large convection 
velocity near the suction side, spreading of the suction leg of the wake segment is much 
faster. When the wake segment reaches the trailing edge of the pressure surface, the wake 
is chopped on the suction side of the leading edge (Fig. 5-16). 

As observed by others, the wake propagation through the passage generates a 
system of the counter-rotating vortices. Analysis of the unsteady flow shows two main 
sources of unsteady pressure on the blade surface. The interaction of counter-rotating 
vortices with the passage flow is the primary source of the unsteady pressure field beyond 
20% of the chord (Fig. 5-17). The second source is the interaction due to pressure gust 
associated with the nozzle wake. A rapid decay of the inlet unsteady pressure field limits 
the influence of the pressure gust to the leading edge region. The interaction of the 
pressure gust with the blade at the leading edge generates a zone of high pressure on the 
suction side (Fig. 5-17). 

Additional simulations were carried out to investigate the contribution of the 
velocity gust and the pressure gust due to the nozzle wake and the rotor blade. In the first 
case, only the velocity defect is present; the second case has only pressure gust. 
Distribution of the instantaneous velocity and pressure fields for these three inlet 
conditions are shown in Fig. 5-20. In the absence of the velocity defect, the pressure gust 
generates a system of counter-rotating vortices with rotation in a direction opposite to 
those induced by velocity gust. Near the leading edge, the variation in the surface pressure 
due to the velocity gust and the pressure gust has a phase angle variation of nearly 1 80°, as 
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shown in Fig. 5-21, at x/C* = 0.71. At this location, the pressure gust has a dominant 
effect. As a result, combined flow has only one peak related to the pressure gust. 

An interesting feature of the pressure gust-blade interaction is that it generates a 
significantly smaller pressure oscillation on the pressure side near the leading edge in 
comparison with the suction side. This can be explained by the shorter interaction time on 
the pressure side. 

In the absence of the pressure gust, pressure oscillations are small within the initial 
15% of the chord. The only significant effect of this interaction is the decreased skin 
friction coefficient on the pressure surface. At moderately off-design flow conditions, this 
phenomenon may lead to an earlier development of the separation bubble due to the 
nozzle-wake blade interaction and may generate significant additional losses. 

Beyond 20% of the chord, the pressure gust has a minimal impact on the unsteady 
pressure field on the suction surface (Fig. 5-22). Downstream of this location, the 
interaction between the nozzle wake induced counter rotating vortices develop a high- 
pressure zone upstream of the nozzle wake center and a low-pressure zone downstream. 
This unsteady pressure system propagates downstream on the suction surface at the 
convection velocity (i.e. velocity at the boundary layer edge near the suction surface). The 
location of the maximum pressure is close to the location of the maximum wake defect. 
Downstream of x/C x = 0.7, the unsteady pressure decays due to a decreased intensity of 


counter-rotating vortices. 
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Interaction between the passing wake and the pressure surface is significantly 
weaker because of a very small variation of the velocity field in the vicinity of the blade 
downstream of x/C x = 0.3. The resulting amplitude of the unsteady pressure on the suction 
surface is less than one-third of those at the pressure surface. 

The measured and predicted unsteady pressure coefficient on the suction surface at 
the design condition is shown in Fig. 5-22. Many of the features were described earlier. 
The measurement does not show the effect of the pressure gust near the leading edge, as 
there is no data available at this location. But the trend from 5% to 20% of the chord 
indicate that the distribution is similar to the predictions. The distribution beyond 20% of 
the chord shows the major influence of the velocity gust and the associated counter- 
rotating vortices. Predictions are in good agreement with the data. Both measurements 
and predictions indicate that unsteady pressures are negligible beyond about 70% of the 
chord. 

5.9 Nozzle Wake Decay Through the Rotor Passage 

The nozzle wake is a source of additional mixing losses in turbine passages. 
Physics of the wake mixing consists of two main components: viscous dissipation and 
inviscid effects (chopping, stretching, distortion, area changes etc.). The wake decay due 
to the viscous dissipation results in losses, while the inviscid effects are the reversible 
process. In some cases, wake ingestion can be used to increase the efficiency of the 
propulsion system (Smith, 1993). 
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Based on experimental and analytical investigations. Hill et al. (1963), suggested a 
wake decay model in a diffusing flow, which includes both the viscous and the inviscid 
effects. Van Zante et al. (1997) developed an inviscid wake model for the two-dimensional 
compressor blade row. This model estimates the inviscid decay of the wake defect as a 
function of the cascade flow parameters. The wake propagation in a turbine is more 
complex in nature when compared to those in a compressor. Due to a strong variation of 
the flow velocity in both the streamwise and the pitchwise directions, the wake is highly 
distorted and can not be considered as straight or slightly bent segments. To analyze the 
wake stretching in a turbine, it should be considered as a set of segments undergoing 
stretching/compression under significantly different local conditions. 

An understanding of the viscous and inviscid contributions to the wake mixing is 
essential for minimization of the mixing losses. An analysis of the wake mixing based on 
the viscous and inviscid prediction is carried out to investigate the physics of the nozzle 
wake decay in a rotor passage. Three following cases are simulated: 

1) Base flow, viscous simulation (includes inviscid effect) 

2) Base flow, inviscid simulation 

3) Modified base flow, viscous simulation; no unsteady pressure variation at inlet. This 
case denoted as ‘no pressure gust’ 

A comparison of the wake decay in the absolute frame of reference for all cases is 
shown in Fig. 5-15. At the trailing edge, the nozzle wake defect based on the inviscid 
prediction is three times larger than that predicted by the viscous simulation. For locations 
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upstream the leading edge in addition to the distribution of the maximum wake defect, 
minimum wake defect distribution is plotted at same locations. These lines correspond to 
the wake propagation along path II, Fig. 5-16, (see discussion below) and denoted la, 2a, 
3a for cases 1,2,3, respectively. 

Viscous and inviscid predictions have a different wake decay rate inside the rotor 
passage, but a similar wake flow pattern exists at all locations. This enables a comparison 
between the inviscid and the viscous wake decay at different locations inside the rotor 
passage, even though the local wake defect is different in viscous and inviscid predictions. 
It is assumed (and verified through the numerical modeling) that at each particular location 
the inviscid rate of the wake decay is not a function of the local wake velocity defect. 
Thus, the inviscid contribution to the wake decay from the viscous solver can be 
compared with those based on the inviscid simulation despite the difference in local 
amplitude of the wake defect. 


Table 5-2 Zone boundaries 


Zone number 

Beginning x/C x 

End x/C x 

1, max AV, corresponds to Path I, Fig. 5-16 

-0.24 

-0.1 

la, min AV, corresponds to Path II, Fig. 5-16 

-0.24 

-0.1 

2 

-0.1 

0.05 

3 

0.05 

0.25 

4 

0.25 

0.6 

5 

0.6 

outlet 

Total 

-0.24 

outlet 

Blade region 

0.0 

1.0 
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Different mechanisms contribute to the wake decay and mixing in a rotor passage. 
There are significant variations in the viscous/inviscid effects at each axial location. It is 
possible to segregate the rotor passage into zones based on the characteristics of the 
nozzle wake decay. The boundaries of the various zones are shown in Tab. 1. The first 
zone extends from the inlet computational plane to x/C x =-10% upstream of the leading 
edge. In this region, the nozzle wake is advanced as a set of parallel segments convected 
by the mean flow. The wake segment, transported along path I (Fig. 5-16) undergoes 
stretching in the streamwise direction, due to a 20% decrease in the streamwise velocity. 
This effect is clearly seen in Fig. 5-15 (line 2), in the form of the increased inviscid wake 
defect. A decrease in the inviscid wake defect beyond x/C x = -0.15 is associated with the 
potential effect, which is described below for the second zone. The wake segment 
traveling along path II (Fig. 5-16) undergoes an opposite process. Flow acceleration (40% 
increase in the local velocity) upstream of the suction surface results in a significant 
decrease in the wake defect (Fig. 5-15). The ratio of the maximum to minimum inviscid 
wake defect at x/C x =0.1 is equal to 1.45. Viscous dissipation has a more profound 
influence on the wake propagation along path I in comparison with the Path n. 

For each zone, the contribution of the viscous dissipation can be calculated as: 
wake decay due to the viscous dissipation 
8 AV d i s = 5 AVvis - 8 AVjnv 

change in the wake defect due to the inviscid effect 8 AV str = 8 AV inv 
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where AV = (V -V)/V is a local wake defect, 6AV - is the change in wake defect between 
the two axial locations, subscripts inv, vis -denote inviscid and viscous solution 
respectively. 

Viscous contributions to the wake mixing, presented in Fig. 5-23 are defined as 
follows: 


and inviscid : 

SAV dis +\SAV inv \ 5AV,„ + |SAVg 

For zone N > 2, the predicted inviscid wake at the inlet to the zone is scaled by the 
predicted viscous wake defect at the current location, i.e., (5 AVj nv )’=5 

A V inv* ( A V vij/ A V inv) . 


g A ^ y/ 

Contributions of the potential effect are calculated as: 5 a V = 

F po ' SAV 


where the subscript visnopot denotes Case 3. A positive sign indicates that the pressure 
gust increases the wake decay. Similar to above, the wake defect scaling is used to 
normalize the wake 8 AVvisnopoi • 

The inviscid mixing is responsible for 30% of the total wake decay along Path I 
(denoted as Zone 1 in Fig. 5-23). For Path II, intense stretching of the wake segment 
increases the inviscid contribution to about 50%. Downstream of x/C x =0. 1, only the 
propagation of the maximum wake defect is considered. 

The essential difference between the first and the second zone (from 10% upstream 
of the leading edge to 5% of the chord from the leading edge) is the increased contribution 
by the pressure gust. Viscous simulation, without the unsteady pressure field at inlet (Case 
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3) predicts a significant increase in the maximum wake defect over the baseline (Case 
3,Fig. 5-15). The comparison, shown in Fig. 5-20 can be used to interpret the influence of 
the moving pressure gust on the nozzle wake decay. This influence can be separated into 
two parts. The first effect is the generation of an additional favorable pressure gradient, 
resulting in enhanced mixing. At the inlet, the wake segment is located in the region of 
high pressure. When the wake is convected downstream, the region of high pressure is 
replaced with a zone of low pressure due to relative nozzle-rotor movement. A numerical 
simulation with “no wake” boundary conditions reveals another contributing factor in 
increased wake decay due to the pressure gust. Even though no time-dependent velocity 
component was presented at inlet, the moving pressure field generates jet-wake structure 
shown schematically in Fig. 5-20. The jet, associated with the moving pressure field, is 
located above the nozzle wake. Thus, the velocity at this side of the wake is increasing, 
contributing to the decay of the wake defect. Counter clockwise rotation of the wake 
segment leads to a compression of the wake segment near the leading edge of the pressure 
side. This counterbalances the decay in the wake defect due to the pressure gust. The 
combination of these two effects results in a smaller net change in the inviscid decay of the 
wake defect. Overall, the contribution of the viscous and inviscid mechanisms is 
approximately equal in this region (Fig. 5-23). 

In the third region (21one 3), which extents from 5% of the chord up to 25% of the 
chord, the location of the maximum nozzle wake defect is moving rapidly from the 
pressure side to the suction side (Fig. 5-16). After the nozzle wake is chopped at the 
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leading edge, two counter-rotating vortices develop in the rotor passage. At x/C x =10%, 
the slip velocity vector (Fig. 5-1) is pointed in the direction of the suction surface. 
Transport of a low momentum fluid from the pressure surface to the suction side is the 
dominant mechanism responsible for the amplification of wake defect in this region. The 
acceleration of the convection velocity in the direction of wake propagation leads to wake 
thickening and a further increase of the wake defect. A viscous dissipation contribution is 
only 60% of inviscid one in this region (Fig. 5-23). 

At x/C x =30% (Fig. 5-15), the wake is strongly bent with two legs extending from 
the suction to the pressure surface. A maximum wake defect is located on the suction 
side, closer to the wake center. The slip velocity is parallel to the streamwise direction. 
Further propagation of the wake segment containing a maximum wake defect can be 
considered through its rotation-free elongation-compression in the streamwise direction. 
Following the velocity field the wake defect is decreasing from x/C x =30% to x/C x = 60%. 
The wake defect increases from x/C x =60% to x/C x =100%, since the streamwise total 
velocity decelerated by 13% (Fig. 5-15). 

In both zones 4 and 5, the inviscid phenomenon is the major mechanism for the 
wake decay (Fig. 5-23) and it is responsible for more than 70% of changes in the wake 
defect. In zones 4 and 5, the wake decay can be estimated using the formula for inviscid 
stretching. The estimated values are as follows: 


Table 5-3 Wake decay, zone 4 and 5 



Predicted 

Inviscid stretching (estimated) 

Zone 4 

0.64 

0.68 

Zone 5 

1.08 

1.20 
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Within the blade passage, the inviscid stretching results in a net increase of the wake 
defect (Fig. 5-15, line 2). In summary, the viscous dissipation is responsible for 60% of the 
wake decay within the blade passage and for 75% of the wake decay from the inlet to the 
trailing edge. 

5.10 Loss Generation Due To Unsteady Flow 

The loss estimate is one of the most challenging tasks for the computational fluid 
dynamists. Numerical simulation based on different numerical solvers may result in an 
appreciable variation in the predicted losses (e.g., Venable et al. 1998). However, the 
prediction of the trend in losses is consistent. Based on a comparison between the 
predicted and measured flow characteristics presented earlier, a certain level of confidence 
is achieved in the predicted loss distribution. However, the absolute level of the predicted 
loss should be verified through a comparison with the measurement. An additional inviscid 
numerical simulation was carried out to assess the presence of the non-physical losses due 
to the numerical errors (artificial dissipation, grid influence, etc.). Grid identical to those in 
the viscous solver was used to ensure consistency. It was found that the time mean (the 
pitchwise mass-averaged) loss coefficient based on an inviscid calculation is less than 0.01 
or less than 4% of the maximum time-mass averaged loss at the outlet based on the 
viscous solver prediction. This is an acceptable level of accuracy. The predicted flow and 
the pressure field can be used to analyze the distribution of losses through the passage as 
well as additional losses due to the presence of the nozzle wake. 



123 


Instantaneous distribution of the relative total pressure, shown in Fig. 5-19, reveals 
that the inlet flow pattern associated with the low total pressure inside the wake is 
followed by a zone of high and low unsteadiness in total pressure along the blade passage. 
Static pressure variation is the major factor contributing to the total pressure oscillation. 
The range of the total pressure variation at the location of the maximum nozzle wake- 
blade interaction (x/C x ~0.3) is significantly higher than the time-average value of At 
this location, the local peak value of C, varies from -1.5 to 0.5 near the suction surface, 
while the time-averaged loss coefficient is equal to 0.25 at the passage outlet (Fig. 5-24). 
The unsteady relative total pressure coefficient reaches its minimum (and correspondingly 
max. amplitude) near x/C x = 0.7. Further downstream, a spot of low total pressure 
propagates downstream without significant change in its value. A zone of high relative 
total pressure undergoes a more rapid decay from x/C x =0.7 to x/C x = 1.5. This is due to a 
more rapid mixing of the nozzle wake. 

The pitchwise averaged loss coefficient, shown in Fig. 5-24, can be used to identify 
sources of increased losses due to the unsteady flow. The loss coefficient for the unsteady 
flow is based on the unsteady inlet total pressure. Thus, it includes only losses in the rotor 

passage. There is a sharp increase in £ upstream of the leading edge; an indication of the 

intense mixing in this region. This contribution accounts for 55% of the increased total 
loss coefficient in comparison with those based on the steady flow prediction. Additional 
losses within the blade passage contribute 34% to the total unsteady losses. The steady 
and unsteady loss distributions are practically parallel from 25 to 50% of the chord. Thus, 
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the profile losses due to the rotor boundary layer are dominant in comparison with the 
mixing losses due to the nozzle wake decay within the passage. The time-averaged 
boundary layer properties beyond x/Cx=0.4 reveal higher values for suction surface 
momentum thickness in comparison with the steady flow prediction. This correlates well 
with the differences in £ based on the steady and the unsteady flow predictions. There is a 
moderate increase ( £ -0.07) in the loss coefficient downstream of the trailing edge. This 
may be attributed to the nozzle wake - rotor wake interaction. However, this level of £ 
increase is to small to be considered as reliable. 

5.11 Unsteady Transitional Flow on the Suction Surface 

An accurate prediction of the unsteady transitional flow is crucial in evaluating 
loses, efficiency, and cooling requirements of turbines. Very few attempts have been made 
to predict the unsteady transitional boundary layer on a turbine blade using the full Navier- 
Stokes solver. 

A major feature of the unsteady flow in the present turbine rotor is the location of 
the maximum unsteadiness in the surface pressure, which is at x/C x =0.3 (s/C x =0.45), (Fig. 
5-22). This is close to the location of the transition onset, thus significantly affecting the 
flow in the transition region. An accurate prediction of the pressure field is crucial to the 
prediction of the unsteady transitional zone. 

Space-time distribution of the unsteady turbulence field, due to the interaction of 
the nozzle wake with blade surface boundary layer, is shown in Fig. 5-25. t/T=0 line 
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corresponds to the flow conditions in Fig. 5-16. Various zones in the unsteady boundary 
layer can be classified as follows (shown in Fig. 5-25 and Fig. 5-26): 

L - laminar region 
TR1 - wake induced transitional strip 
TU 1 - wake induced turbulent strip 
TR3 - transition between wakes 

TR2- wake turbulence induced transition between wakes 

TU2 - turbulent flow between wakes 

TR4 - transitional strip behind passing wake 

The unsteady nature of the transitional process makes it difficult to identify the 
beginning and the end of the transition. For a steady flow simulation, the identification of 
the inception and the end of the transition can be based on the local turbulence field and 
verified through an analysis of other flow characteristics (i.e., skin friction coefficient, 
shape factor, etc.). It is found that this approach cannot be used in an unsteady flow due 
to a phase lag between different flow variables and properties in the transitional region. 
Hence, the location of the onset and the end of the transition is based only on the 
predicted turbulence field. The position of the transition inception is calculated as the 
position where F((i t m ) reaches its maximum value along the path on a space-time diagram 
(Fig. 5-25). F((i t m ) is the first derivative of |i, m along the wake path at the boundary layer 
edge and (i t m is the maximum value of the eddy viscosity inside the blade boundary layer 


at any given space-time position. 
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The location of the maximum value of k at the end of the transitional zone was 
chosen as the location of the end transition. This criteria correlates well with the predicted 
value of the local- maximum skin friction coefficient. 

The variation of the predicted unsteady skin friction coefficient on the suction 
surface is shown in Fig. 5-26. As indicated earlier, from x/C x =0 to x/C x = 0.25 the pressure 
gust is the primary source of the unsteady pressure in the blade (Fig. 5-22). 

Beyond x/C x = 0.2, the wake defect and the counter-rotating vortex system 
generate a combination of high-pressure and low-pressure zones, which are convected in 
the streamwise direction. In the laminar region, the surface pressure variation is the 
primary source of the boundary layer disturbance. The observed variation of the skin 
friction coefficient along Path 1-2, la-2a as well as 3-4, 3a-4a, shown in Fig. 5-26 can be 
explained on the basis of the unsteady pressure gradient (favorable and adverse) at these 
locations (Fig. 5-22) 

Near x/C x =0.3 (s/C x =0.45) an interaction of the blade boundary layer with the 
nozzle wake increases the velocity at the boundary layer edge (Fig. 5-16). At this location 
(Fig. 5-22 and Fig. 5-25), the Reynolds number, based on the momentum thickness, (Ree,) 
is 12% higher over the time-averaged value, reaching a value of Ree=230. The boundary 
layer development along this path (going through point 1 in Fig. 5-26) undergoes wake 
induced transition at s/C x =0.68. Transition within the wake induced transitional strip is 
characterized by a maximum increase in the turbulence kinetic energy as well as the 
maximum value of the unsteady skin friction coefficient. After the boundary layer becomes 
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fully turbulent, the nozzle wake disturbance leads to a wake induced turbulent strip (TUI). 
The transformation of the kinetic energy of the velocity fluctuation into the turbulence 
kinetic energy generates an increased level of k across the boundary layer in comparison 
with the zone of the turbulent flow between wakes (TU2). Declining pressure fluctuations 
and wake defect decay are the reasons for the vanishing shear stress fluctuations in the 
wake induced strip beyond 80% of the chord. 

A laminar flow in the region following the wake has a moderately decreased level 
of Re e =200 at 30% of the chord (Fig. 5-26, Point 2). The boundary layer development 
behind the nozzle wake is characterized by a significantly delayed transition inception, 
which occurs at s/C x =0.98. Flow remains transitional up to the trailing edge (Zone TR4). 
Unlike the wake induced transition strip or the transition region between wakes, the 
turbulent kinetic energy does not undergo a rapid increase in the transition region, with 
further relaxation to the level of those in the developed turbulent boundary layer. The 
turbulent kinetic energy increases monotonically; raising from a low laminar value to those 
in the end of the transitional zone. An analysis of the boundary layer development in this 
region leads to the conclusion that the diffusion of the turbulence kinetic energy between 
the ffeestream and the turbulent zones upstream and downstream are the primary 
mechanisms responsible for the amplification in k. 

The zone in between the passing wake is not a subject .to an additional pressure 
gradient or other disturbance associated with the wake-blade interaction. Nevertheless, 
there is an earlier transition inception near s/C x =0.70 and t/T=0.5. A comparison of the 
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turbulent kinetic energy (Fig. 5-25), and the blade-to-blade turbulence field (Fig. 5-18) 
shows an elevated level of turbulence in the boundary layer due to the maximum diffusion 
of the nozzle wake turbulence at this location. In the absence of the velocity disturbance, 
this high turbulence level (about 5 % increase in Tu, ) is the primary source of earlier 
transition in Zone TR2. A numerical simulation with no turbulence fluctuation in the 
nozzle wake confirmes this hypothesis. If there is no elevated level of turbulence in the 
nozzle wake. Zone TR2 moves upstream and closer to Zone TR1, merging with Zone 
TR3 (transition between wakes). 

The boundary layer in the Zone TR2 has a lower increase in the momentum 
thickness and has a lower level of eddy viscosity. Thus, losses in this region are lower in 
comparison to those at the wake induced transition/turbulent strips, but a little higher than 
that in the region of pure transition between wakes (Zone TR3). 

The flow in the Penn State turbine has a high level of freestream turbulence. For 
this type of flow, the transition occurs through the bypass mechanism. The transition 
inception occurs at the location where an outer disturbance causes generation and growth 
of turbulent spots. The length of the transitional region depends on the spot spreading 
angle and their growth. These characteristics are functions of a number of parameters; 
local momentum thickness, pressure gradient, Mach number, etc. At zero pressure 
gradient, the leading edge of the turbulent spot propagates at O. 88 W 5 , while the trailing 
edge has the propagation velocity equal to O. 5 W 5 . The low Re k-e model, originally 
developed for the fully turbulent flow, models the transitional process through the 
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diffusion of the ffeestream turbulence and a local balance between the turbulence 
production and the dissipation rate, which depend on the distribution of the velocity in the 
boundary layer. The unsteady flow results in the unsteady transition on a blade surface. 
Disturbances, caused by the passing wake, lead to an earlier generation of the turbulent 
spot, which may form a turbulent strip. Turbulent spot growth and propagation 
mechanisms are essentially the same as in the case of the steady transition. The 
propagation velocities of the leading and the trailing edge of the wake induced turbulent 
strips are equal to those of the turbulent spot (i.e., 0.88Wg and O. 5 W 5 . correspondingly). 
This is confirmed by the measurement (e.g., Halstead et al., 1995). Development of the 
calmed region, which is located behind the wake induced transitional/turbulent strips, is 
another crucial element in the unsteady transitional process. The trailing edge of the 
calmed region propagates at 0.3Wg. This region has laminar characteristics, but has an 
elevated level of shear stresses. This can be considered as the relaxation zone between the 
turbulent and the laminar zones. Higher resistance to the separation and low boundary 
layer losses associated with a calmed region are beneficial. 

The low Re k-e model lacks the physics associated with the turbulence spot 
development and can model the transition process only in an “global fashion”. In the 
numerical prediction, the unsteady velocity fluctuation caused by the wake-blade 
interaction is the primary source responsible for the earlier transition. The maximum value 
of the average turbulence kinetic energy and the maximum value of the k fluctuation 
across the boundary layer occur at the same location. After the transition inception, two 
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major mechanisms are responsible for the further development of the predicted transitional 
strip - the convection of the transitional zone at the local propagation speed, which is 
smaller than the wake propagation speed, Wg at the boundary edge. The second 
mechanism is the advection of the velocity disturbance with variable velocity across the 
boundary layer. It was found that the propagation speed of the transitional strip TR1, 
defined as a path along maximum k in Fig. 5-25, is equal to 0.66Wg. This is higher than 
the propagation velocity equal to 0.55Wg at the location of the maximum turbulent kinetic 
energy across the boundary layer. Investigation of the time-space distribution of the 
velocity and turbulence across the boundary layers reveal that this increase is due to a 
greater speed of the velocity disturbance propagation; the second contributor to the 
development of the unsteady transitional zone. From Fig. 5-25 and Fig. 5-26 the leading 
edge velocity of the transitional strip is estimated to be 0.9-l.Wg , while the trailing edge 
velocity is about 0.5-0.55Wg. As stated above, the transition in Zone TR2 is due to an 
elevated level of ffeestream turbulence, which is converted at the local propagation speed. 
Therefore, the propagation speed of the transition strip is equal to the local convection 
velocity inside the boundary layer (0.55Wg). 

Zone B (Fig. 5-25 and Fig. 5-26) has several features of the calmed region. A low 
level of the turbulent kinetic energy and eddy viscosity in this region is similar to those 
observed in a laminar flow. At the same time, it has an elevated level of shear stresses in 
comparison with the shear stress in the laminar zone. Phase lag between the maximum 
amplitude of k and the maximum amplitude of Cfat the transition inception (about 5% of 
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t/T, Fig. 5-25 and Fig. 5-26) results in higher shear stresses in the zone above the wake 
induced transitional strip. This phase lag diminishes further downstream and is equal to 
zero at the end of the transition. There is a significantly delayed transition onset in Zone B. 
A low level of the momentum thickness in combination with low eddy viscosity results in 
significantly smaller boundary layer losses. Nevertheless, this region lacks another crucial 
feature of the calmed region. A predicted shear stress profile does not have “relaxing” 
distribution of Cf at the side of the calmed region. Experimental data indicate that the 
trailing edge of the calmed region is propagating at O.3W5. In the prediction, the line of a 
zero momentum thickness variation, which can be considered as a boundary of the region 
B, has a propagation velocity equal to 47% of Wg. 

5.12 Rotor Wake Development 

Numerical modeling accurately predicts the decay of the rotor wake defect (Fig. 
5-27). Both the experimental data and the numerical simulation show a rapid decay of the 
rotor wake velocity defect. Analysis of the instantaneous velocity field (Fig. 5-28), as well 
as the wake space-time distribution at x/C x =1.32 (Fig. 5-16), may be used to understand 
the rotor wake development. 

The development of the rotor wake is influenced by its interaction with segments 
of the nozzle wake. The nozzle wake propagation and its distortion in the rotor passage, 
discussed earlier, result in a zone of decelerated flow, which is located near the suction 
surface. A similar zone, with significantly smaller amplitude, propagates along the pressure 
surface (Fig. 5-16). Due to the differences in the convection velocity, the phase lag 
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between these two spots is t/T == 0.35 at the trailing edge. Downstream of the trailing 
edge, the interaction between the rotor wake and these two nozzle wake segments is the 
primary source of rotor wake unsteadiness with a variation in the suction side velocity 
playing a dominant role. Analysis of the pitchwise distribution of the total velocity reveals 
the presence of the rotor wake (Fig. 5-28) and the suction side, and pressure side nozzle 
wake segments within about 5% of the chord downstream of the trailing edge (Fig. 5-28). 
Downstream of x/C x = 1.05%, the suction side wake segment merges with the rotor 
wake. The phase lag between the suction side segment and the pressure side segment of 
the nozzle wake gives a rise to fluctuation in the rotor wake at the double nozzle wake 
frequency, but with a significantly smaller amplitude. 

A comparison of the predicted rotor wake development with the LDV 
measurements shows that the rotor wake development is simulated correctly (Fig. 5-27). 
The merging of the suction side nozzle wake segment with the rotor wake and the 
presence of the pressure side wake segment far downstream is simulated correctly (Fig. 5- 
28). 

Even though a good correlation is achieved for the rotor wake velocity defect, the 
current solver underpredicts the wake semi-width. There is a discrepancy between the 
predicted and the measured wake turbulence levels in the wake region. Measured 
turbulence field shows a very complex three-dimensional anisotropic nature of the rotor 
wake turbulence. Hence, it can not be adequately modeled using the k-e model. 
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5.13 Off-Design Flow 

The numerical simulation at the off-design condition (1100 Rpm) has been carried 
out to derive an additional insight on the unsteady flow physics. Reduced frequency for 
this condition is Cl = 7.32. 

Only the unsteady blade surface pressure distribution has been measured for this 
operating condition. Due to a lack of the experimental data on the inlet boundary 
conditions, the wake profile at the design condition has been utilized for the numerical 
simulation. Velocity and pressure distribution, shown in Fig. 5-29, indicates a significant 
change in the overall flow pattern at the off-design condition. Due to an increased reduced 
frequency, there are more nozzle wake segments in the rotor passage in comparison with 
the design case. Another significant factor is the development of a strong separation 
bubble on the pressure surface, which extends from x/Cx = 5% to x/Cx =65%. In the zone 
of the maximum thickness (x/Cx = 0.35) the separation zone occupies nearly 25% of the 
rotor passage. There are two factors leading to differences between the design and off- 
design conditions. The unsteady pressure field on the suction surface does not undergo 
significant adjustment beyond those caused by an increased reduced frequency (Fig. 5-29). 

The interaction between the nozzle wake and the separation bubble results in 
amplified unsteady pressure on the pressure surface. Predicted unsteady velocity (Fig. 5- 
29) indicates that the nozzle wake generates a vortex pattern in the separation zone. 
Additional pressure fluctuation associated with the vortex structure is the source the 
increased variation in the surface pressure. Based on the space-time distribution, it is 
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possible to identify three major regions on the pressure surface. From the leading edge to 
x/Cx = 0.20, vortex generation and its downstream propagation leads to a zone of 
high/low pressure regions propagating at local convection speed. A comparison of the 
experimental and measured pressure distribution at x/Cx = 0.14 (Fig. 5-11) shows that the 
numerical prediction was able to capture this phenomenon. For x/Cx > 0.2, the 
propagation of the vortex structure does not directly affect the surface pressure. The 
phase change is negligible. This can be seen in the measured pressure distribution (Fig. 5- 
10). A comparison of the measured and predicted pressure at x/Cx =0.3 (Fig. 5-11) 
reveals a phase lag between the measured and the predicted data. This is an indication that 
the numerical simulation overpredicts the extent of the separation bubble. This may be due 
to several factors: 1) possible three-dimensionality of the flow in this region; 2) deficiency 
of the k-e turbulence model for the case of a strong separation bubble. Downstream of the 
separation bubble (x/C x = 0.65), the distribution of the unsteady pressure is similar to 
those at the design condition. 

An overall comparison between the measured and the predicted pressure field 
shows that the numerical simulation predicts the unsteady pressure field fairly well. 

5.14 Effect of Inlet Turbulence 

Additional simulation studies have been carried out with different freestream 
turbulence intensities at inlet. The design condition had a 7% turbulence intensity in the 
free stream at the inlet. The additional simulation includes a freestream turbulence 
intensity of 2% and 15% (which is closer to that existing in aircraft engines). The major 
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feature is the decay of the nozzle wake. The amplitude of unsteady pressure distribution 
shown in Fig. 5-30 reveals that the freestream turbulence intensity plays a major role. The 
wake decay is slower and unsteadiness is higher at the low turbulence intensity and 
decreases significantly as the turbulence intensity is increased to 15%. This is confirmed 
by the nozzle wake defect plotted in Fig. 5-3 1 . 
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Fig. 5-1 LDV and Kulite measurement locations in rotor passage (velocity 

triangle is based on design) 
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Fig. 5-2 Variation of rotor inlet static and stagnation pressure coefficients, 
velocity and turbulence quantities at x/c x =-0.244, design conditions (case A) 



Fig. 5-3 Decomposition of the revolution periodic data into blade periodic (average 
passage) and blade aperiodic 


137 



Fig. 5-4 Averaged blade pressure and unsteady pressure envelope at design conditions, ( 
midspan ) 
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Fig. 5-6 Blade-to-blade distribution of relative velocity 
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Fig. 5-8 Measured space-time distribution of blade unsteady pressure, case A (design ) 



Fig. 5-9 Measured and predicted harmonics of blade unsteady pressure suction surface 
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1100 RPM 



Fig. 5-10 Measured space-time surface pressure distribution at midspan, case C 
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Fig. 5-1 1 Measured and predicted time history of unsteady pressure, case C (off design) 
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Fig. 5-13 Nozzle wake defect at the nozzle-LDV position N1 
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Fig. 5-15 Wake defect decay 
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Fig. 5-17 Unsteady pressure field inside the rotor passage 
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Fig. 5-18 Transport of the nozzle wake turbulence in the rotor passage 







x/Cx=0.071 


boundary conditions 

combined 

— only pressure gust 




x/Cx=0.571 



x/Cx=0.85 






Prediction 



Fig. 5-22 Unsteady pressure coefficient, suction surface 
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Fig. 5-25 Maximum turbulence kinetic energy in the boundary layer, suction surface 
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Fig. 5-26 Unsteady skin friction coefficient, suction surface 
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Fig. 5-27 Rotor wake decay as a function of streamwise 
distance distance 


rotor wake at x/Cx=1 .32 W/W 1 



Fig. 5-28 Space-time distribution of the total velocity at x/C x = 132% 
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Fig. 5-29 Unsteady pressure and velocity fields, case C 
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Fig. 5-30 Variation of the maximum blade unsteady pressure at different inlet turbulence 
levels 






Chapter 6 


NUMERICAL SIMULATION OF THE TRANSITION OVER LAMINAR 

SEPARATION BUBBLE 


The flow in a low pressure turbine stage of the aircraft engine is complicated by 
the variation of the flow conditions. A significant change in the ambient condition results 
in notable variations in the stage Reynolds number. Modified boundary layer development 
affects losses, efficiency, and heat transfer characteristics of the stage. At cruise condition, 
the flow Reynolds number may be less than half of the value of the take-off condition. This 
may result in a separated flow and efficiency degradation. The development of a reliable 
prediction technique is very important in the design of efficient machines and may lead to 
an improved efficiency and weight/thrust characteristics. The transition in a low-pressure 
turbine may occur in either bypass form, similar to those observed in high pressure turbine, 
investigated in the previous chapter, or through the development of a separation bubble, 
depending on the Reynolds number between the take-off and the cruise condition. 

Considerable effort has been spent in the investigation of the ability of different 
turbulence models to predict various types of transitional flows. Nevertheless, very few 
investigators were focused on turbomachinery flows with the transition over a laminar 
separation bubble, especially at a high level of freestream turbulence (Michelassi et al., 
1997, transonic turbine , Huang and Xiong., 1998, low pressure turbine). 
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The objective of the research presented in this chapter is to gain a detailed 
understanding of the unsteady transitional flows in low-pressure turbines, with emphasis 
on separation-induced steady transition. The test case chosen for this study is the 
simulation of the separation and the transition of the flow over the suction surface of low- 
pressure turbine cascade blade investigated experimentally by Qiu and Simon (1997). The 
influence of the free stream turbulence and pressure gradient is investigated. An existing 
Navier-Stokes unsteady flow solver is used. Three low-Reynolds number forms of two- 
equation turbulence models have been incorporated and tested for accuracy. In order to 
overcome the over-prediction of the turbulence kinetic energy and the dissipation rate near 
the leading edge, several modifications of the production terms have been incorporated in 
the code as well as the Algebraic Reynolds Stress Model. 

6.1 Description of the Test Case 

The experimental data in a simulated LP turbine cascade have been used to assess 
the ability of the numerical solver. A schematic of the facility is shown in Fig. 6-1 (Qiu 
and Simon, 1997). The cascade flow was simulated using a channel with a convex and 
concave wall profiled as suction and pressure surfaces of a turbine blade. A flow suction 
device was utilized to simulate periodic flow near the leading edge. Experiments were 
carried out with the inlet flow velocity ranging from 3 to 12.5 m/s, which corresponds to 
Re number from 50, 000 to 200, 000. A number of turbulence generators were utilized to 
generate the flow with 0.5, 2.5 and 10% inlet turbulence intensities. Boundary layer 
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characteristics were measured using a hot-wire probe. Coordinates of measurement 
locations are given in Table 6-1. According to Qui and Simon (1997), the uncertainty in 
the mean velocity is 3.6%, and the fluctuating velocity is 4%. 

Table 6-1 Location of the experimental data points on the suction blade surface 


N 

x/C x 

N 

x/C x 

PI 

0 . 

P8 

0.6889 

P2 

0.0398 

P9 

0.7457 

P3 

0.2111 

P10 

0.8173 

P4 

0.3778 

Pll 

0.8593 

P5 

0.4667 

P12 

0.9111 

P6 

0.5506 

P13 

0.9728 

P7 

0.6247 




6.2 Numerical Procedure 


Numerical simulation shows that for the LP turbine flow, the utilization of the 
explicit ARSM does not modify the solution outside the boundary layer, but may cause a 
stability problem. To avoid this instability and minimize the CPU time, the ARSM is used 
only up to twice the boundary layer thickness from the blade surface. 

Experimental data were acquired assuming two-dimensionality of the flow. 
However, turbulent and especially transitional fields are essentially three-dimensional. To 
address this feature, the numerical simulation based on the ARSM/ k-e approach is carried 
out using a three-dimensional solver. The symmetry boundary condition in the spanwise 
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direction forces the mean flow to be two-dimensional. Nevertheless, the turbulence field 
contains the spanwise component of the Reynolds stresses. 

Experimental data were acquired in a “channel” type configuration, while an 
original “cascade” configuration was chosen for the numerical simulation. The utilization 
of the bleeding device to model the cascade flow may lead to a discrepancy in the flow 
angle in the vicinity of the leading edge. To verify the potential effect of this discrepancy, a 
numerical simulation of the cascade flow at different inlet flow angles have been carried 
out. A comparison of the predicted and the measured blade pressure distribution indicates 
that the best prediction has been achieved at the design flow angle (Fig. 6-2). 

Another potential source of the discrepancy between the “channel” and “cascade” 
configurations is the flow near the trailing edge. An extension wall employed in the 
experiment provides a smooth development of the boundary layer beyond the point of the 
virtual trailing edge. In contrast, the flow over a real trailing edge is characterized by a 
sudden change in the flow angle, local variation of the pressure, etc. The numerical 
experiment shows that, in the case of the attached flow, utilization of either approach 
leads to a practically identical solution except in the region 5 % upstream of the leading 
edge. However, the flow prediction based on “cascade” configuration results in significant 
instability as soon as the reattachment point moves downstream of the trailing edge. In 
order to reduce the effect of this phenomenon, but keep the “cascade” approach, the blade 
has been extended using an extension wall with zero thickness for cases with a separation 


bubble. 
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The numerical investigation requires a verification to ensure the grid independence. 
A number of grids (121 x 71, 141 x 91, and 241 x 181; or 121x71x3 in the case of 
ARSM/k-e model etc.) have been utilized to study the grid dependency on the solution. 
The maximum distance between the surface and the first grid point varies between y+=0.8 
for the course grid to y+=0.12 for the finest grid. At x/c x =0.65, coarse grid has 20 grid 
points within the boundary layer and 8 grid points within the laminar sublayer. The fine 
grid has 35 and 15 points correspondingly. Numerical predictions based on coarse and fine 
grid are very close to each other. In some cases, the fine grid solution was not stable in the 
transition region over a separation bubble. The difference between the solutions based on 
fine and coarse grid is minimal for the converged solution. All the reported simulation data 
in this paper is based on 141 x 91 grid. Additional discussion on the code verification is 
presented later. 

The numerical simulation of the transitional flow in turbomachinery cascades 
reveals the deficiency of the standard k-e model in predicting the flow with a high free 
stream turbulence. Non-physical increase in the turbulence intensity near the stagnation 
point may “contaminate” the boundary layer turbulence and trigger an earlier transition. 
An elevated level of the turbulence at the mid-passage, in the zone of maximum flow 
acceleration, leads to a higher level of ffeestream turbulence (ATu= 3-4%) at the boundary 
layer edge. The modifications of the k-e model described in Chapter 3 is used to predict 
turbulence filed near the leading edge and at midpassage accurately. After the 


incorporation of the production term modification, the inlet distribution of the dissipation 
rate is set to provide the measured kinetic energy along the boundary layer edge. 


6.3 PREDICTION USING k-e MODEL 


6.3.1 Case Re=200000, Tu=10% 

The distribution of the surface pressure predicted by different turbulence models is 
compared with the experimental data in Fig. 6-2. Since the flow is fully attached at this 
Reynolds number, the pressure increases monotonically along the rear part of the suction 
surface. There is no difference between the blade pressure distribution predicted by 
various turbulence models. Flow with Re = 200,000 and Tu=10% is fully attached on the 
suction surface (Fig. 6-3). The flow with Re=50,000 and Tu=10% has transition over 
laminar separation bubble as shown in Fig. 6-4. 

In the laminar part of the boundary layer (experimental locations P2-P7, Table 6- 
1), the predicted velocity field exactly match the measured values. For brevity, this 
comparison is not shown. Velocity and turbulence intensity profiles at locations P8-P13 
are shown in Fig. 6-5. The beginning and the end of the transition, as well as the 
separation location predicted by various turbulence models, are compared with the data in 


Table 6-2 
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Table 6-2 Inception and length of the transition, separation and reattachment points, 
Re=200000, Tu=10%. 



Experiment 

Prediction 
FLB model 

Prediction 
CH model 

Prediction 
LB model 

Transition inception, x/C x 

71% 

-65% 

57-59% 

61% 

End of transition, x/C x 

81% 

82% 

94% 

80-82% 

Separation, x/C x 


attached 

75% 

attached 

Reattachment, X/C x 

— 


no reattachment 



1 Flow visualization indicates the presence of a very small separation bubble at x/C x =0.7 


In the laminar boundary layer, the prediction based on the LB model is identical to those 
based on the FLB model. The turbulence intensity profile shows that the CH model 
predicts an amplification of the turbulent kinetic energy in the laminar boundary layer. 
While the eddy viscosity predicted by the CH model is four-five times higher than the eddy 
viscosity derived from the simulation based on the FLB or LB models, it is still small. The 
ratio of eddy viscosity to laminar viscosity is about 4 and this part of the boundary layer 
has laminar features. The primary source of this increase is the lack of turbulent 
dissipation to compensate for an increase in the turbulent kinetic energy near the leading 
edge, rather than generation inside the boundary layer. The increased level of turbulence 
intensity predicted by the CH model contributes to an earlier transition. 

The simulation based on the CH model predicts the transition at about 58% of the 
chord. The best agreement between the measured and the predicted velocity profile is 
achieved in simulations based on FLB and LB models. At locations P8 and P 1 1 the 
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predicted velocity profiles are identical to the experimental data. Between P8 and PI 1 the 
measured profile is less steep in comparison with the numerical solution. The transition 
inception is predicted at about 4% of the chord upstream of the measured location. The 
numerical simulation predicts a shorter transition and hence the flow field is near the end 
of the transition at the location P10. The turbulence profile is very close to fully 
developed’ profiles. As a consequence of this, the predicted velocity profiles at P9-P10 are 
closer to turbulent profile when compared to the data. 

The main drawback of the CH model is the existence of a very thin separation zone 
near the 75% of the chord. The CH model, based on y+, is known for its poor 
performance in separated flows. The CH model overpredicts the turbulence intensity in the 
turbulent boundary layer by 50-70%. 

6.3.2 Case Re=50000. Tu=10% 

Reduced Re number results in the development of a medium size separation bubble 
(Fig. 6-4). The separation zone is characterized by the presence of a flat zone in the 
pressure distribution at x/Cx~0.75 on the suction surface of the blade (Fig. 6-2b). A 
solution based on the FLB model correctly predicts this trend, which shows an earlier 
return to an adverse pressure gradient. As confirmed by the velocity profiles, this is due to 
a smaller separation bubble and earlier reattachment of the flow. In spite of the presence of 
the separation bubble in the simulation based on FLB and CH models, the region of 
constant pressure is not clearly predicted by these two models. 
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Velocity profiles are plotted in Fig. 6-6. Experimental data indicates that the flow 
separates at about 74% of the chord and the transition to turbulence occurs over a laminar 
separation bubble at about 8 1 % of the chord. All three models predict a laminar separation 
at the same location; x/Cx=71% of the chord. The length of the laminar region inside the 
separation bubble varies for different models. The CH model gives the transition 
immediately after the inception of the separation at 71% of the chord. A simulation based 
on the FLB model predicts the transition inception at 78% of the chord, with distance 
between the separation point and inception of transition equal to 7%, which is close to the 
measured value. 

The transition to turbulence over the separation bubble is characterized by the 
inception of the transition in the shear layer with further penetration through the 
separation zone. The maximum turbulence intensity is located farther from the wall in 
comparison with high Re cases (Fig. 6-6). Contrary to the measured data, the numerical 
simulation predicts strong backward flow inside the separation flow, (point P10). An 
additional turbulence production due to the higher shear stresses in the zone of separated 
flow increases the turbulence intensities in the separation bubble. The predicted turbulence 
profile has a smoother distribution, and its maximum is located closer to the wall in 
comparison with the experimental data. In the case of FLB and LB models, an increased 
level of turbulence in the separation zone leads to a smaller thickness of the separation 
bubble and an earlier reattachment. A comparison between the numerical simulation and 


the data is given in Table 6-3. 
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Table 6-3 Inception and length of the transition, separation and reattachment points, 
Re=50000, Tu=10%. 



Experiment 

FLB model 

CH model 

LB model 

Transition inception, x/C x 

76% 

78% 

72% 

75% 

End of transition, x/C x 

>pl3 

97.28% 

84% 

85% 

86% 

Separation, x/C x 

70% 

71% 

72% 

72% 

Reattachment, X/C x 

91% 

93% 

no 

reattachment 

no 

reattachment 


6.3.3 Case Re=50000. Tu=2.5 % 

This case is the most difficult to compute. A low level of turbulence at a low Re 
number leads to an inherently unsteady flow with an unsteady separation bubble and 
transitional zone. Even though the results presented in this paper is based on the steady 
solution, the analysis of the convergence and unsteady flow simulation indicate the need 
for an implementation of the time accurate simulation to achieve better resolution of the 
flow physics. No results for the LB model is presented, because attempts to stabilize the 
solution using increased artificial dissipation resulted in total damping of the separated 
flow. For the CH and FLB models, the flow has moderate fluctuation in the size and 
extent of the separation bubble. The data presented is calculated as an average of these 
fluctuations. This variation affects the rear part of the separation region and does not 
influence the location of the separation point and the transition inception point. 

A comparison between the predicted and the measured surface pressure 
distribution is shown in Fig. 6-2c. Both models underpredict the extent of the separation 
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bubble, which results in a shorter zone of constant pressure. A comparison of the 
measured and the predicted velocity and turbulence intensity profiles (Fig. 6-7) shows 
similar trends to the cases described earlier. The numerical solver overpredicts the 
turbulence intensity in the transition zone; while for the flow with high freestream 
turbulence level, the maximum turbulence intensity is underpredicted (similar to 
Re=200000 case with Tu=2.5% and Tu=10%). There is no peak in turbulence fluctuations 
above the separation bubble in the transition region. The prediction based on the FLB 
model has a smaller separation bubble thickness and an earlier reattachment. The size of 
the separation bubble is equal to the experimental data at point P9, about 1/2 of the 
experimental value at point P10 and 1/3 at point PI 1. For the FLB model, the variation in 
the thickness of the separation bubble was 50%. A low level of the freestream turbulence 
delays reattachment from the 85% to 99% of the chord. A comparison between the 
numerical simulation and the data is given in Table 6-4. 


Table 6-4 Inception and length of the transition, separation and reattachment points, 
Re=50000, Tu=2.5%. 



Experiment 

Prediction 
FLB model 

Prediction 
CH model 

Prediction 
LB model 

Transition inception, x/C x 

81% 

-80% 

-78% 

69% 

End of transition, x/C x 

96% 

86% 

88% 

97% 

Separation, x/C x 

69% 

71% 

69% 

72% 

Reattachment, X/C x 

no 

reattachment 

99% 

no 

reattachment 

90% 
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6.4 Prediction Using Hybrid ARSM/K-e Model 

Experimental data for the transition over a laminar separation bubble (e.g., Wang 
and Hatman, 1998) shows a strong redistribution of the turbulent kinetic energy between 
components in the transition zone. Peaks of u’ and v* are equal to each other, while in the 
attached turbulent boundary layer v’ is equal to 40% of the u’ component. The k-£ model 
is unable to capture this redistribution as well as the overall anisotropy of the turbulence 
field associated with the transition process. A numerical simulation based on the hybrid k- 
e/ARSM has been carried out to investigate the ability of this model to improve the 
prediction of the transition flow over the LP turbine blading. Results of the current 
research as well as previously reported simulations (e.g., Abid et al., 1995) indicate that a 
numerical solution strongly depends on the k-£ model used. A comparison of the 
prediction based on hybrid models with different low Re k-£ (CH, FLB, LB) led to the 
conclusion that transition inception is controlled by the k-£ model and is close to those 
predicted by a corresponding k-£ model. Therefore, the FLB model has been chosen as the 
model with the best results based on the previous computations. 

Results of the numerical simulation based on hybrid k-£/ARSM are identical to the 
prediction based on k-£ approach. High Re case is characterized by the transition in an 
attached boundary. Maximum shear stresses are located close to the blade surface in the 
region where viscous residual is calculated using k-£ approach. ARSM part is used only 
for the outer layer. Therefore, the influence of the hybrid approach is minimal. In contrast, 
the shear layer above the separation zone is located in the zone where viscous residual is 
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calculated using ARSM. Thus, the hybrid approach has a more profound effect. For the 
high Re case (Re=50 000, Tu=10%), this effect is minimal and can be seen only in the 
turbulence field. The utilization of a hybrid model moves the peak of the fluctuation 
velocity further from the wall, and closer to the measured location. There is only a minor 
change in the predicted velocity field. A comparison between the predicted and the 
experimental data for the case with Re=50, 000 and Tu=2.5% (Fig. 6-8) reveals an 
improvement in both predicted velocity and turbulence. 

Current numerical simulations have been carried out without the pressure strain 
terms. As a result, the w’ component is equal to the v’ component. An analysis of the 
turbulent case indicates that the streamwise component has about 50% of the total 
turbulent kinetic energy, while v’ and w’ have 25% each. No significant change in the 
balance between the different components is found in the transition region (no more than 
5% variation). Hybrid turbulence model overpredicts maximum amplitude of the 
fluctuation velocity similar to k-e model. However, redistribution of the turbulence energy 
between turbulence components plays a major role in improving the velocity prediction . 

6.5 Prediction Using k-e Model in Conjunction with the Transition Model 

A numerical simulation of the transitional flows based on the turbulence model 
generally does not provide an adequate level of accuracy and robustness. Incorporation of 
transition models is a potential way to improve the transition prediction. Transition models 
use an empirical or a semi-empirical correlation to calculate the inception and end of 
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transition, as well as the intermittency distribution in the transition region zone. A number 
of models were developed to calculate the inception and end of transition in an attached 
flow. The most common approach is the calculation of the transition inception using an 
empirical correlation and the calculation of the intermittency distribution using the 
approach suggested by Dhawan and Narasimha (1958), in conjunction with the correlation 
for the non-dimensional spots breakdown parameter (e.g., Gostelow and Walker, 1991; 
Mayle, 1991). In the current research, a model by Abu-Ghannam and Show (1980) is 
utilized for the transition prediction in an attached flow: 

Re e „ = 163 + exp |V(A) - - ( ^~ — 


TO) = 


J6.9 1 + 1 2.75 • A + 63 .64 • A 2 • • ifk < 0 
[6.91 + 2.48 • A - 12.27 • A 2 • ■ ifk > 0 


k = S£_ 

where: ^ e 

transition is completed when: 

Re e = 2Re ftr 

Intermittency distribution is based on Dhawan and Narasimha (1958) 

Transition models for separated flows correlate the distance between the 
separation and the transition inception. In the current paper the model due to Davis et al. 
(1985) is used: 

Re„ r - Re^ = 25 • 10 3 log 10 [coth(17.327u e )] 
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As in the case of the attached flow, an intermittency distribution is based on Dhawan and 
Narasimha (1958) formula. Even though this relation was suggested for the attached flow, 
a comparison presented in Qui and Simon (1997) indicates that it can be applied to the 
separated flow. 

Utilization of the transition model provides a more reliable prediction in 
comparison with the “pure” turbulence model. Nevertheless, the transition model has a 
number of weak points. Empirical correlations are based on data that has a significant 
spread. The accuracy of the transition model deteriorates as flow parameters deviate from 
those used for the derivation of the model. A comparison between the predicted 
distribution and the measured data shows that the transition model used predicts an earlier 
transition. Therefore, this discrepancy may be attributed to the inaccuracy of the transition 
model. In addition to the models described earlier, calculations with the transition model 
based on the measured intermittency factor have been carried out to investigate an effect 
of the “ideal” transition model. For the attached flow case (Re=200, 000, Tu=10%), there 
is no significant difference in the predicted velocity and turbulence fields because the 
predicted transition inception based on Abu-Ghannam and Shaw correlation is located 
upstream of the “natural” transition inception predicted by a “pure” turbulence model. For 
the separated flow, the transition model used is practically identical with the experimental 
data for the case with Re=50, 000 and Tu=10%. However, as described in Qui and Simon 
(1997), for other cases with a high freestream turbulence, the correlation was not perfect. 
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Different approaches can be used to incorporate an intermittency distribution. 
Methods based on a separate solution for the laminar and turbulent parts show potential 
for the improved flow prediction (e.g., Steelant and Dick, 1996). Nevertheless, for 
engineering applications it is more advantageous to have a single solver throughout the 
flow field. This is especially preferred from the point of view of the model extension to 
three-dimensional flows in turbomachines. Two methods are tried for the incorporation of 
the intermittency factor into the existing code. In the first method, an additional damping 
function F (y) is utilized for the calculation of the eddy viscosity (this case of the transition 
model utilization is denoted as Var.l): 

M,=— c„/„F(r) 

£ , where F(y)=y 

This approach implicitly assumes that the eddy viscosity based on local scales is 
‘turbulent’ in nature. If the distribution of the turbulent kinetic energy and the turbulent 
dissipation rate are transitional in nature, this approach may lead to an underprediction of 
the local eddy viscosity due to the ‘double’ damping. Assuming that the eddy viscosity in 
the transition region can be calculated correctly using the same expression for (it as for the 
turbulent part, the intermittency distribution can be utilized only through the modification 
of the turbulent equations. An additional damping function based on y is applied only to 
the calculation of the production term (this case is denoted as Var. 2): 

P = P * F(y) ^ w h ere p i s the turbulence production 
Mean flow is affected implicitly through the eddy viscosity. 
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In the transition zone, the intermittency distribution has complex two-dimensional 
distribution. Most of the transition models vary intermittency only in streamwise direction. 
Based on experimental distribution, the calculations have been carried out with the 
variation of the y in both streamwise and crossflow directions. Overall six combinations of 
transition model incorporations have been analyzed; there are two ways to incorporate 
intermittency distribution into solver (Var.l and Var. 2) and three different distributions of 
the intermittency factor: 

1) Step distribution: y=0, for x<Xtr and y =1, for x>Xtr 

2) One dimensional :y(x) = Max y (y(x,y)) or :y(x) base on transition model 

3) Two-dimensional : y^y(x,y) 

All cases are calculated using the FLB turbulence model, which gave the best prediction 
among the CH, LB, and FLB models. 

In attached flow transition, the implementation of the transition model does not 
have any significant effect on the velocity and turbulence distribution. No significant 
influence of the method of the implementation of the transition model (Var. 1, Var. 2, or 
type of y) is found. As stated above, the utilization of Abu-Ghannam and Shaw correlation 
predicts an earlier transition in comparison with both the experiment and the prediction 
based on “pure” turbulence model. The current approach may only postpone the transition 
inception. Therefore, the prediction based on the “pure” turbulence model and the 
transition model produce practically identical flow fields. Simulation with the experimental 
distribution of the intermittency factor improves the prediction the turbulence kinetic 
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energy distribution in the vicinity of the transition inception and an earlier end of the 
transition. Nevertheless, there is no improvement in the velocity distribution at location 
P10. 

In contrast to the high Re cases, the way the transition model is incorporated and 
the type of intermittency distribution used has a profound effect on the prediction of low 
Re flow (Fig. 6-9). The incorporation of the transition model with the direct effect on the 
eddy viscosity (Var. 1) resulted in the development of the larger separation bubble in 
comparison with the experimental data (Fig. 6-9). The separation zone extends beyond the 
location of the trailing edge. Utilization of the two-dimensional distribution of the 
intermittency factor led to a further increase in the separation bubble size. The flow 
prediction based on the application of the intermittency distribution to the calculation of 
the production term (Var.2) leads to the prediction of a much smaller separation bubble 
and reattachment near the trailing edge. The predicted height and extent of the separation 
zone are closer to the measured values, which is a consequence of the delayed inception of 
transition. However, an overall deviation of the predicted velocity profile from the 
experimental data is greater in comparison with the “pure” turbulence model for all cases 
except the “step” transition model. This is due to the double damping of the eddy viscosity 
in a transition zone. Even though y distribution indicates that the transition zone should 
extend beyond the trailing edge, all but one (one-dimensional model, Var. 1) has the end 
of the transition upstream of the trailing edge. A numerical prediction based on the one 
dimensional distribution and “step” distribution in conjunction with Var. 2 gave the most 
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accurate prediction of the separation bubble size and location, even though it does not 
improve the turbulence intensity distribution in comparison with the simulation without the 
transition model. 

6.6 Effect of Artificial Dissipation on the Transition Prediction 

The stability consideration requires an analysis of the differential approximation of 
the original PDE, which is not possible. Numerical experiments are carried out to study 
the influence of the artificial dissipation on the transition prediction. The major objective 
of the current validation is to establish the limits of its influence. 

A number of realizations of the artificial dissipation terms have been analyzed. The 
original version of the code employed a hybrid second/fourth order artificial term with a 
switch based on the local turbulence field. Velocity scaling and eigenvector scaling are 
incorporated to keep the artificial dissipation at a minimum level in the boundary layer. 
Nevertheless, an analysis of the turbulence kinetic energy in the transition zone (Fig. 6- 
10b) indicates that the level of the artificial dissipation reaches 50% of the source term (Pk 
- e) for the base case. The base case has k2ke=001, which is about twice the minimum level 
required to avoid odd-even decoupling. To minimize the level of the artificial dissipation, 
the artificial dissipation term was modified to include only the streamwise variation of k. 
The k balance based on this modification is shown in Fig. 6- 10c. All calculations 
presented in this paper are based on this modified approach, even though it does not affect 
the prediction beyond a small zone of high gradients in the transition region. 
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The effect of utilizing only the fourth-order artificial dissipation term in turbulence 
transport equations has been also investigated (for mean flow equations the artificial 
dissipation is always based on fourth-order terms). This approach does not alter the result 
of the analysis presented below (beyond absolute values of the artificial dissipation 
coefficient). However, the employment of only the fourth-order artificial dissipation leads 
to a significantly increased sensitivity of the code to the turbulence field development near 
the leading edge. A moderate flow disturbance generates a significant increase in the 
turbulence kinetic energy, which decreases rapidly downstream. Numerical modeling 
shows that this increase can not be explained as a transition with relaminarization further 
downstream, because it may be reproduced at any location within the first 30% of the 
chord by placing the source of potential disturbance (e.g., locally skewed grid). 

The predicted location of the separation inception, beginning and end of transition 
and reattachment point as a function of the artificial dissipation is shown in Fig. 6-11. 
Both the insignificant as well as excessive levels of artificial dissipation result in an earlier 
transition. The values of k4 and k 2 ke vary from the level below the stability limit to a level 
at which the artificial dissipation causes a significant non-physical diffusion. 

It should be noted, based on the previous experience with the solver, that the 
recommended variation of the k 2 ke was 0.01 — 0.02. Within this range, the variation of the 
predicted and measured location of the transition inception is within 2.5% of the chord. 
An earlier transition inception results in a smaller separation bubble (Fig. 6-1 1 and Fig. 6- 
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12). However, the variation of the artificial dissipation in the mean flow equation does not 
significantly influence the predicted mean flow profile (Fig. 6-13). 

The primary source of the early transition in the case of a small k 2 ke is a slight 
numerical instability of the scheme. For k 2 ke ^ 0.075, a moderate odd-even decoupling 
generates a premature transformation from the laminar to the turbulent boundary layer. An 
increase in the artificial dissipation also results in an earlier transition inception. It is 
possible to identify zones with a different behavior of the scheme. For simulations with 
k 2ke <0.02, the variation of the artificial dissipation term affects only the transition 
inception, but the transition length is essentially constant. This fact indicates that, within 
this range, the artificial dissipation acts as a destabilizing factor. A comparison of the 
streamwise distribution of the turbulent kinetic energy based on differing values of k 2ke 
shows that slope of k is constant; i.e., the transition zone is shifted upstream without 
diffusion of the k field. Therefore, in this zone, the artificial dissipation is similar to the 
physical disturbances (freestream turbulence, noise etc.). For k2 ke >0.02, the artificial 
dissipation leads to both an earlier transition and an increased transition length. This is the 
consequence of the streamwise/stream diffusion of the turbulent kinetic energy. 

The influence of the artificial dissipation in the mean flow equation on the 
predicted velocity and turbulent fields is presented in Fig. 6-13. In contrast to the k 2 ke 
variation, the variation of k 4 does not affect the accuracy of the mean flow, except at very 
high levels of artificial dissipation. 
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6.7 Concluding remarks 

Comparison between the prediction and the experimental data formally shows 
more accurate prediction of the transition inception for cases with the transition over a 
separation bubble. However, the development of the attached, transition and turbulent 
boundary layers is not very sensitive to even relatively high error(~4 % of the chord) in the 
predicted transition inception. In contrast, the development of the separation bubble is 
strongly affected by the small variation in the transition prediction. Moderate delay in 
transition inception results in significant enlargement of the separation bubble and 
correspondingly notable increase in profile losses. To be reliable, the numerical solver 
should provide the prediction of the transition inception with an accuracy equal to 1-2% of 
the cord. Further research is needed to achieve this goal. Numerical investigation, 
presented in this Chapter, shows very strong effect of solver characteristic on the 
transition development for the flow with the transition over a separation bubble. 
Therefore, further development should be concentrated not only on the transition 
model/turbulence model development but also on the model coupling with the numerical 
solver to control the influence of the scheme characteristic on the transition development. 





Fig. 6-2 Surface pressure distribution 
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Fig. 6-9 Velocity and turbulence intensity, Tu=10%, Re=50000, effect of the transition 
model utilization a) Id model Varl, b) Id model Var2, c) “step” model, d) 2d model 
Varl, e) 2d model Var2 
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Fig. 6-10 Turbulent kinetic energy balance, Tu=10%, Re=50000 
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Fig. 6-11 Effect of artificial dissipation on predicted location of transition and separation 
inception, , Tu=10%, Re=50000 
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Fig. 6-12 Influence of the artificial dissipation; Tu=10%, Re=50000 
k4=0.01 a) k 2ke =0.001 b) k 2ke =0.004 c) k 2ke =0.02 d) k 2ke =0.04 



Fig. 6-13 Influence of the artificial dissipation, Tu=10%, Re=50000 
k 2ke =0.01 a) k4=0.001 b) k4=0.02 c) k4=0.04 



Chapter 7 


NUMERICAL SIMULATION OF LEADING EDGE FILM COOLING 

Numerical simulation of two dimensional steady and unsteady transitional flows 
presented in previous chapters provides a foundation for accurate modeling of 
turbomachine aero thermodynamics. Turbine film cooling flow and secondary flows, 
including tip gap flow are two major problems that have profound effect on the 
characteristic of turbomachinery performance. The ability of the numerical solver to 
predict these phenomena should be established. Numerical simulation of the three- 
dimensional leading edge film cooling flow, discussed in this chapter, has been performed 
to accomplish this objective and gain a better understanding of the vortex structure due to 
the cooling jet-main flow interaction. The heat transfer and the flow phenomena associated 
with the leading edge film cooling are very complex. A thin boundary layer near the 
stagnation point, large pressure gradient, and the presence of the curvature effect make 
the numerical modeling including the grid generation extremely complicated. This is 
further complicated by the compound angle injection. Very few attempts have been made 
to simulate such flows. The research presented in this chapter is limited only to the steady 
state approximation of the leading edge film cooling due to the large CPU resources 
required for adequate temporal resolution of the unsteady phenomenon. Multiblock 
version of the solver, described in Chapter 2, has been utilized for the numerical modeling. 



191 


This has been found to be crucial for the accurate prediction of the flow in film cooling 
multidomain configurations. 

7.1 Flat Plate Film Cooling 

The film cooling on a flat plate with an injection from a row of holes, measured by 
Pietrzyk et al., (1990) was chosen as a case to test the code and establish a numerical 
procedure for the numerical simulation with cooling flow injection. A similar configuration 
was computed by Leylek and Zerkle (1993) (incompressible code) and Fabian (1995) 
(ADPAC developed by Hall et al., (1994)). The inlet flow velocity was 20 m/s. Blowing 

ratio M = — -SiSS l . = 2 . Ratio of the hole length to the hole diameter 1/d was 3.5. Inlet 
(pU) 

mean flow 

turbulence intensity was T u =0.5%. The computational domain was divided into 3 blocks: 
plenum, hole, and main part. Grid (Fig. 7-1) was generated using Genie+ grid generator in 
conjunction with an algebraic grid generator for the grid point distribution inside the holes. 
The grid size was 81x31x31 for the main block, 31x17x17 for the hole, 21x33x33 for the 
plenum. To ensure a stable and converging calculation, the numerical simulation was 
carried out through the number of steps. During the first step a flow initialization was 
carried out in separate blocks. Two dimensional calculation was performed to generate 
initial flow distribution in the main block. Multidimensional flow in a cooling hole was 
calculated using constant outlet static pressure. Hole inlet pressure was set to provide the 
required massflow. The next step consist of a coupled calculation of the flow in the 
plenum and the hole, preserving the prescribed blowing ratio. During the third step, the 
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flow in the mainstream was calculated with a fixed jet velocity from the previous step to 
avoid the numerical instability and reverse flow from the main block into the hole. Finally, 
a coupled simulation of all three blocks was carried out to obtain a solution. 

The result of the numerical prediction was compared with the experimental data 
and the numerical prediction based on the ADPAC code. The predicted and the measured 
distribution of the adiabatic effectiveness along the jet centerline is shown in the Fig. 7-2. 
Due to the high blowing ratio, the jet separates from the surface upstream of the cooling 
hole. This leads to a sharp decrease in the adiabatic effectiveness at s/d=3. The comparison 
of other flow characteristics (secondary vectors, crossflow temperature distribution) 
indicates a good correlation between the prediction and available data. 

7.2 Leading edge cooling at compound angle 

The configuration and experimental data by Cruse et al. (1997) are utilized for the 
numerical modeling of the leading edge film cooling aimed at validation and improved 
understanding of the flow and thermal physics. The schematic of the model is shown Fig. 
7-3. The symmetry of the flow near the leading edge was simulated using the suction slot 
near the stagnation point. The cooled air (T= 166 K) was injected through two rows of 
holes. The first row was located along the stagnation line and the angle between the hole 
axis and the surface was 20°. The second row of holes was inclined at 25° to the 
upstream direction. The angle between the hole axis and the surface was 20°, which is the 
same as the corresponding angle for the first row. There were nine holes in each row, 
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with the hole spacing equal to z/d=7.67. The length of the hole was equal to l/d=12. The 
upstream plate and the suction were adjusted to obtain a correct position of the stagnation 
streamline. The surface temperature distribution was measured using an infrared camera. 
Thermocouples were employed to measure the crossflow temperature distribution. 

7.2.1 Numerical modeling 

Based on the results of the preliminary analysis, the multiblock grid, presented in 
Fig. 7-4, was generated for a numerical simulation. It consists of 6 blocks: a plenum, 2 
holes, a mainblock, an inlet, and an outlet block. The total number of grid points is 
286199. A summary of the flow condition used is presented in Table 7-1. 


Table 7-1 Flow conditions 


Freestream Velocity (m/s), Uo 

10 

Mass Flow/ Hole (g/s) 

0.725 

Freestream Turbulence (%), Tu 

0.5, 20. 

Pressure (atm),p 

1 

Freestream Temperature, (C),T 0 

27.5 

Plate Conductivity (W/mK) 

0.025 

Average Mass Flux Ratio, M 

2.0 

Surface Roughness (m) 

<25 

Jet/ Freestream density ratio, D 

1.8 

Hole Edge Radius (mm) 

<0.1 


The excessive grid skewness may result in a nonphysical solution due to the 
increased level of the artificial dissipation. For the current topology, the angle between the 
hole wall and the cooling wall is equal to 20°. In the case of the structured grid used in 
the current simulation, the presence of the sharp angle requires an optimal redistribution of 
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the local skewness near the hole-main block interface. It was found that the accuracy 
could be improved if some non-orthogonality is allowed near the hole walls as opposed to 
the case with the orthogonal grid near the hole wall and a higher concentration of skewed 
cells farther from the wall. The boundary conditions used in the numerical modeling are as 
follows (letters correspond to block faces shown in Fig. 7-4): 


Table 7-2 Boundary conditions 


ai 

Inlet (P\ T\ flow angle) 

a 2 

inlet (p,U, flow angle) 

b 

Periodic conditions 

c 


d, 

no-slip conditions 

d2 

wall function 

e 

outlet, constant static pressure 


Interaction between the coolant jet and the mainflow results in a very complex 
flow field. The three-dimensional separation zone downstream of the cooling jet has an 
adverse effect on the overall heat transfer. This complexity requires the implementation of 
the low-Re number turbulence model instead of the wall function. While the use of the 
wall function can significantly reduce the computational time, it lacks the near wall physics 
required for an accurate resolution. To ensure an accurate prediction of the boundary layer 
flow in the main block, the first grid point was located at y +s =l .2. 

The static pressure at the plenum inlet was set to provide the required massflow 
through the injection holes. The current experimental configuration has long holes, 
l/d=12. It has been found that after the developed flow regime is achieved near the hole 
entrance, the plenum flow does not affect the flow in the mainblock. Therefore, the 













195 


plenum flow field is frozen during the final convergence to minimize CPU utilization. The 
pressure oscillation during the convergence process resulted in the development of the 
reverse flow into the hole from the mainblock. To avoid this and ensure a stable 
development of the jet, the pressure distribution for the hole is fixed during the initial part 
of the computation. After a fully developed jet is obtained in the mainblock, this limitation 
is relaxed. 

In the first attempt, a low-Re turbulence model was also employed for the flow in 
the hole. Conformal (one-to-one) interblock interface requires that stretched grid lines 
from the near wall region inside the cooling hole should be extended throughout the main 
block. Thus, the quality of the mainblock grid is significantly reduced due to the presence 
of the zones with very stretched (in normal to the cooling surface direction) and skewed 
grid cells. Consequently, the numerical solution has suffered from an excessive numerical 
dissipation; and the predicted mixing of the coolant jet with mainflow has been strongly 
overestimated. The incorporation of the wall function for the boundary condition inside 
holes has enabled the generation of the grid with a more uniform distribution of the grid 
points on crossflow planes. The mainblock grid also has smoother characteristics (as a 
result of one-to-one interblock interface). Despite a minor sacrifice in the accuracy of the 
flow prediction in the coolant hole due to the utilization of the wall function, the accuracy 
of the mainblock flow prediction has improved significantly. 
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7.2.2 Comparison with the experiment 

A comparison between the predicted and the measured surface adiabatic 
effectiveness is shown in Fig. 7-5. Smooth distribution of the adiabatic effectiveness in the 
vicinity of the stagnation line indicates the uniform spreading of the coolant flow from the 
bottom row of jets. In the numerical simulation, a re-circulation zone exists between the 
jet and the wall and; as a result, the coolant jet separates from the wall. Interaction with 
the main flow causes the reattachment of the jet at point A (Fig. 7-5a). The extended 
region of the cooled surface near point B (Fig. 7-5b) in the experimental data is an 
indication of the similar behavior (reattachment) of the jet. More intensive mixing and 
diffusion of the bottom jet in the experiment result a in a low surface temperature 
between upper holes (Fig. 7-5,point C). In the numerical prediction, the bottom coolant 
jet encircles the root of the top jet and is partially entrained by it. This contributes to the 
predicted low adiabatic effectiveness at point C, Fig. 7-5. 

The lateral distribution (line plot) of the adiabatic effectiveness due to the presence 
of the upper cooling hole is plotted in Fig. 7-6 at several streamwise locations. The 
numerical simulation accurately predicts the amplitude and the position of the maximum 
cooling at s/d = 4.86. There are two main discrepancies between the predicted and the 
measured data; the predicted temperature decrease on the left side of the jet is less steeper 
than that in the experiment and, the predicted cooling effectiveness is lower than the 
measured values between z/d = 4.5 and z/d = 7. The difference between the predicted 
location of the maximum influence of the bottom jet is one of the factors contributing to 
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this discrepancy. As indicated above, in the numerical prediction the bottom jet reattaches 
far to the left (A, Fig. 7-5a) as compared to the experiment (B, Fig. 7-5b). As a result, the 
rest of the coolant fluid from the bottom jet augments cooling on the left side of the top jet 
rather than in the region between z/d= 4 and 7. At s/d= 9.88 (Fig. 7-6c), the minor 
discrepancy between the predicted and the measured maximum adiabatic effectiveness 
indicates an undertuming of the predicted jet trajectory. The discrepancy at lateral sides of 
the jet, observed at s/d=4.86, practically vanishes due to the mixing process. 

The overall effectiveness of the film cooling process is represented by laterally 
average adiabatic effectiveness. A very good agreement between the measured and the 
predicted laterally averaged adiabatic effectiveness is achieved (Fig. 7-7). Even though a 
good agreement with the experimental (thermal) data was achieved on the surface, the 
numerical simulation underpredicts the temperature diffusion in the core of the jet. This 
characteristic can be found in all numerical simulations of film cooling. Under some 
conditions, this may lead to a more favorable prediction of the surface temperature 
distribution in comparison with the real configuration. A comparison of the predicted 
spanwise distribution of normalized temperature with measured values (Fig. 7-8) indicates 
that the normalized temperature in the core of the jet is 0.25 (s/d= 1.24) to 0.4 (s/d= 9.88) 
higher than the measured values. The isotropic nature of the turbulence model used is 
probably the major source of this problem. The turbulence field associated with the jet 
mainstream interaction is strongly anisotropic. k-£ turbulence models may not be suitable 
for the prediction of such flows. Another factor is the choice of the inlet turbulence 
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dissipation rate. An inlet level of e has less influence on the wall bounded flow, but may 
have a stronger influence on the free shear layer flows. Further evidence of this is 
presented later. 

7.2.3 Discussion of aerothermal flow physics 

A thorough understanding of the physics of the jet-mainstream interaction can 
contribute to an improved understanding of the primary sources of film cooling losses as 
well as factors contributing to the cooling effectiveness. The development of vortices due 
to the jet-mainstream interaction and its effect on mixing has a profound influence on the 
overall development of the film cooling flow, including the distribution of the coolant fluid 
and aerodynamic losses. The importance of tracking the vortex structure is emphasized by 
Sgarzi and Leboeuf (1997), and Walters and Leylek (1997). Sgarzi and Leboeuf (1997) 
suggested a classification of major vortices associated with the film cooling on a flat plate. 
In the case of the flat plate, there are five vortices: “kidney shaped” counter-rotating 
vortices in the core of the jet; “horse shoe” vortices due to the sudden deceleration of the 
boundary layer upstream the leading edge of the jet; “half wake” vortex pair in the zone of 
low pressure downstream of the jet.; “half wall” vortex pair induced by the “kidney” 
vortex; and “lip” vortex, due to the ffeestream jet- leading edge interaction. Visualization 
by Bario and Beral (1996) confirms the existence of these vortices. Walters and Leylek 
(1997) analyzed the source of different vortices as well as their influence on the jet 
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injected over a flat plate. Based on this analysis, suggestions were made to improve the 
film cooling. 

Leading edge film cooling through a compound angle injection results in a much 
more complex aerothermal field in comparison with the flat plate injection. Complexities 
arise due to the difference in the direction of the jet velocity and mainstream velocity, the 
three-dimensional turning of coolant jet, and the presence of strong streamwise pressure 
gradient. 

The numerical solution is analyzed to examine the vortex structure, including the 
origin of various vortices and their contribution to the flow and thermal field development. 
A schematic representation of major vortices associated with the upper row of holes is 
shown in Fig. 7-9 - Fig. 7-12. Streamlines are plotted like ribbons; their twist corresponds 
to the magnitude of the local vorticity. The laminar boundary is very thin upstream of the 
jet. The amplitude of the “horse shoe” vortex (denoted as Q], Fig. 7-9) is very weak. 
Streamlines downstream of the injection hole clearly indicate the presence of the “half 
wake” vortex (denoted as Q 2 . Fig. 7-9). The vorticity in the core of the jet induces a 
counter-rotating pair of vortices (denoted as Q 3 ) which can be seen in Fig. 7- 13a. The 
development of the “kidney” vortex in the core of a flat plate cooling jet was shown to 
have a first-order effect on the film cooling effectiveness. In the case of the leading edge 
film cooling at compound angle, a similar, but more complex, structure exists in the film 
cooling jet. It is possible to identify four main vortices. During the initial part of the jet 
path they are associated with the streamlines from the hole orifice. These streamlines are 



200 


shown in Fig. 7-9 - Fig. 7-12. The major part of the cooling fluid is encountered in vortex 
£Zib,.Fig. 7-11. This vortex rapidly dissipates (Fig. 7-13a, Fig. 7-14a, Fig. 7-15a). 
Downstream of s/d= 1 1 it changes its sign. The cooling fluid associated with this vortex 
comes from the core of the jet. The part of the coolant fluid located closer to the right 
lateral side of the hole is organized into the vortex £Xu. The characteristic feature of this 
vortex is its high strength. At a distance far from the jet (Fig. 7- 13c), it has a structure 
similar to those observed in a flat plate film cooling. The interaction between the jet and 
the mainflow also produces two counterclockwise rotating vortices Fig. 7-11 and 
£2 5b , Fig. 7-10. They have their origin at the upstream side of the hole (Qsb) and left lateral 
side of the hole (Q 5a ). The behavior of these vortices is similar to the clockwise rotating 
counterparts. The vortex Q 5b decays rapidly, while vortex Q 5a retains its strength far 
downstream (Fig. 7-13). 

An understanding of the origin of the vortices is the basis for the improved jet- 
mainflow interaction. Walters and Leylek (1997) identified three potential sources of the 
“kidney vortex,” (the boundary layer at the lateral side of the coolant hole, an axial 
secondary vorticity, and the lateral crossflow shear layer). The first one has a major 
influence on the development of the ‘kidney’ vortex. The analysis of the flow near the 
upper hole in the current simulation leads to a similar conclusion. Vorticity due to the 
boundary layer in the hole at sides A and B (Fig. 7-12) is the primary source of Qsa and 
Q 4 a vortices, as well as £l* b - Vorticity due to the boundary layer at the upstream side of 
the hole (C, Fig. 7-12) contributes to the vorticity fi 5 b- Crossflow shear layers 


201 


predominantly affect the development of vortices, Qsb and Qib- To analyze the 
contribution of different vortices to the overall interaction, the projection of the 
normalized vorticity in the main stream direction, a normalized temperature, and the loss 
coefficient for three different streamwise locations are plotted in Fig. 7-13, Fig. 7-14 and 
Fig. 7-15. The plane at s/d =4.2 is located near the upper hole. At this location, the jet 
undergoes a rapid turning from the crossflow direction to the mainflow as well as bending 
around the z axis. All four of the core vortices defined earlier are clearly seen in Fig. 7- 
15a. The distribution of the normalized total pressure has two major spots. The region in 
the core of the jet has a positive value of £ . The stagnation region behind the jet has a low • 
level of the total pressure (Fig. 7- 15b). The presence of the low pressure zone, which 
extents 2.5d downstream of the trailing edge of the jet (Fig. 7-16) is the primary source of 
the strong crossflow, which can be seen upstream of the stagnation region (Fig. 7-17 and 
Fig. 7-18). The crossflow is compromised of the hot fluid and creates a zone of high 
temperature under the coolant jet. The combined effect of £2s», and CLu vortices and 
crossflow lifts and moves the zone of a low total pressure to the right (Fig. 7- 14b). The 
pressure gradient also results in larger turning of the upper jet in comparison with the 
bottom jet. The experimental results indicate a higher skeweness of the thermal field than 
that observed in the numerical prediction at the left side of the jet at this location. It may 
be an indication that the numerical simulation undepredicts the strength of the Qsb vortex. 
Intensive mixing and dissipation due to the presence of Qsb and il»b vortices decrease the 
value of the total pressure in the core of the jet. At s/d=6, the jet has practically completed 
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its turning in the streamwise direction. Further downstream, the development of the jet is 
similar to that of the flat plate jet, except for the presence of the streamwise pressure 
gradient. The core of the jet reaches the shortest distance to the wall at about s/d=7; 
further downstream it lifts slowly. 

One of the negative effects of counter-rotating vortices is an entrainment of the hot 
air under the cooling jet. Fig. 7-9 - Fig. 7-12 show that £2 5a and vortices entrain hot 
fluid under the jet, pushing cold core fluid farther from the wall, hence decreasing the 
cooling efficiency. An analysis of the streamlines shows that £Xi a and 0 5a are a mixture of 
coolant and mainstream flows. Upstream of s/d~5, £2s a consists of hot air with the rest of 
the coolant fluid pushed closer to the jet center. £X,b and £2 5b have a minimal entrainment 
of the surrounding hot fluid and predominantly consist of a cooling fluid with an increased 
temperature due to the diffusion. 

The contribution of various vortices is summarized below: 

1) The vortex flu has a strong impact on both the thermal efficiency (entrainment of 
the hot flow) and aerodynamic losses. 

2) The vortex Q 5a has a similar effect. Its contribution to the losses is smaller in 
comparison with those of £I» a . The presence of the stagnation region downstream of 
the jet is the primary source of the losses at the left side of the jet. 

3) Vortices £ljb and Q 5 b have a much less adverse effect due to the rapid decay, 
caused by the strong dissipation near the root of the jet and a minor lifting of the jet 


core downstream. 
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4) The three-dimensional flow field arising from the second (upper hole) jet is shown 
in Fig. 7-17 and Fig. 7-18. The presence of the separation zone underneath the jet 
as well as the highly skewed free shear layer above the separation zone is clearly 
visible. The flow angle changes dramatically from the wall to the free stream. Such 
large and differing flow angles of the jet and the mainstream result in large velocity 
gradients and intense turbulent mixing. The jet-wake type of profile can be clearly 
seen in Fig. 7-17 and Fig. 7-18. The mixing of the two jets is completed in about 
six jet diameters. Large crossflow (spanwise flow) develops at the bottom of the 
jet, resulting in the thermal field observed (Fig. 7-17 and Fig. 7-18). These 
crossflows tend to lift the separation zone and move the hot spots away from the 
wall. 


7.2.4 Mach number effect 

All measurements of the film cooling configurations are performed at a very low 
Mach number, usually around 0.03-0.06. The Mach number in a real configuration is 
usually much higher. To analyze the effect of the Mach number on the leading edge film 
cooling, a numerical simulation of the configuration with M=2 and stagnation point 
located at s/d=0 is carried out with an inlet Mach number of 0.3. All other variables used 
for the test case are held constant. Maximum Mach number at the top of the cylindrical 
leading edge reaches 0.9. Due to the rapid acceleration of the flow, static temperature of 
the flow above the surface is much smellier than that at a low Mach number. To make a 
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proper comparison between the high and low Mach number cases, the wall temperature in 

the case of the high speed flow is adjusted by a change in the kinetic energy at the edge of 

2 _ 2 

the boundary layer due to the local acceleration AT = M/ "*— 1 — — UlowMac h . a comparison of 

2C p 

the corrected adiabatic effectiveness between the experiment and the numerical prediction 
at a low Mach number (Fig. 7-19) shows a mild increase in the jet diffusion at a high speed 
flow. Two major factors contribute to this behavior of the jet; a higher pressure gradient 
(Fig. 7-20) and the decreased diffusion of the jet due to the compressibility. Analyses of 
the flow near the stagnation line indicate that there is a moderate decrease in the 
separation zone under the bottom hole jet. Hence, the mixing process is increased due to 
the smaller lifting of the jet. The flow pattern in the vicinity of the top row of coolant 
holes does not indicate modifications of an aerothermal field due to the higher pressure 
gradient resulting from the increased Mach number. Consequently, the adiabatic 
effectiveness is not changed. Further downstream, the jet in the high-speed flow 
undergoes a rapid spreading. This is due to the influence of the negative pressure 
gradient. 


7.2.5 Influence of the high inlet turbulence 

The level of inlet turbulence in a real turbine is much higher than that used in 
laboratory experiments. To understand the influence of the high ffeestream turbulence 
intensity on the leading edge film cooling, a numerical simulation of the flow with 10% 
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inlet turbulence was carried out. This is close to the turbulence intensity encounted in real 
engines. A comparison of the laterally averaged adiabatic effectiveness for high and low 
freestream turbulence intensities (Fig. 7-21) indicates a decrease in the cooling effect near 
the bottom row of coolant jets. Cooling effect downstream of the top row of jets is 
modified only slightly. 

The overall vortex structure resulting from the injection from upper holes is similar 
to that observed for a low turbulence case (Fig. 7-22a and Fig. 7-22b). An increase in the 
turbulence intensity leads to a more intense mixing and more rapid decay of the vortex 
system and results in a significantly low temperature of the jet core (Fig. 7-23a and Fig. 7- 
23b). The position and the amplitude of vortices are also modified due to a more intense 
turbulent diffusion. The most affected vortex is £2^ (Fig. 7-9 - Fig. 7-12 and Fig. 7-13a). 
Its amplitude is significantly reduced and is located closer to the wall at Tu=10% (Fig. 7- 
22a and Fig. 7-22b). This results in the smaller lifting of the £2^ vortex, containing the 
core of the coolant jet. The reduced intensity of £Xu vortex also decreases an entrainment 
of the hot ambient fluid under the coolant jet, therefore improving the cooling 
effectiveness. A distribution of the surface adiabatic effectiveness (Fig. 7-24) shows more 
intense spreading in the crossflow direction. This favorable effect of the increased 
turbulence (low entrainment and higher spreading) compensates for the intensified jet-core 
temperature diffusion. Therefore, the laterally averaged adiabatic effectiveness is modified 
only moderately (Fig. 7-5a and Fig. 7-24a). 
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Increased turbulence intensity significantly affects the flow structure in the vicinity 
of the coolant holes near the stagnation surface. The separation zone under the bottom jet 
shrinks by a factor of two in both the crossflow and in the normal to the wall direction 
(Fig. 7-25 and Fig. 7-26). The location of the jet core moves closer to the wall at a higher 
turbulence intensity. Another major modification is a smaller spreading of the jet in the 
streamwise direction. This is due to the lower intensity of the “kidney”-type vortex. At the 
low turbulence intensity level, the “kidney”-type vortex spreads to the lateral side of the 
jet before turning smoothly in the streamwise direction. At high inlet turbulence, coolant 
jet streamlines go parallel to the stagnation line until the jet turns suddenly in the 
streamwise direction; and without reattachment, finally mixes with the mainstream. The 
concentration of the jet streamlines along the stagnation line compensates only slightly for 
the increased turbulence diffusion. 


7.2,6 Influence of the inlet turbulence length scale. 

The development of the leading edge film cooling is controlled by the vortex 
structure located outside the boundary layer. In boundary layer flows, the freestream 
turbulence has a major influence on the flow development through its influence on the 
turbulence intensity at the edge of the boundary layer. Due to the intense mixing and an 
entrainment of the ambient fluid into the coolant jet core, the freestream turbulence 
dissipation rate and the corresponding length scale play a significant role in the leading 
edge film cooling effects. 
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To analyze the effect of the freestream turbulence length scale, a numerical 
simulation was carried out with an inlet turbulence intensity T u = 10% and a higher length 
scale in comparison with the previous case. For the current calculation, the ratio of the 
eddy viscosity and the molecular viscosity is 120, while for the previous simulation (with 
T„ = 10%) it is about 30. Smaller dissipation rate leads to only a minor increase in the 
turbulence intensity (about 1.5% near the bottom coolant hole row). The comparison of 
the lateral averaged adiabatic effectiveness at two differing length scales is shown in Fig. 
7-27. A comparison of the surface adiabatic cooling effectiveness (Fig. 7-24) reveals a 
more intense increase in the surface temperature at the lower turbulence dissipation rate. 
The increased eddy viscosity moderately reduces the laterally averaged adiabatic 
effectiveness downstream of the upper hole. Examination of the aerothermal field indicates 
only a minor decrease in the intensity of the £X} a vortex (Fig. 7-22b and Fig. 7-22c). 
Therefore, there is no mechanism to compensate for the higher temperature diffusion 
similar to the previous case. 

As indicated by Cruse et al. (1997), the experimental investigation gave a 
contradictory conclusion on the effect of the increased turbulence intensity. The 
modification of the mixing process due to the change in the length scale can be one of the 
factors contributing to these discrepancies. Elevated turbulence dissipation rate affects the 
jet core and correspondingly the vortex strength £lib- A numerical simulation with a low 
freestream turbulence intensity overpredicts the temperature of the coolant jet core (Fig. 
7-8). This simulation was carried out with a high freestream turbulence length scale. The 
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observed trend due to the variation in the length scale (Fig. 7-23b and Fig. 7-23c) 
indicates that the inlet length scale employed for the baseline computation contributes to 
the observed discrepancy between the measured and the predicted jet core temperature. 

7.3 Concluding remarks. 

The numerical investigation indicates that the CFD analysis can be successfully 
employed for the prediction, simulation, and parametric study of the complex flows 
associated with the leading edge film cooling. However, great care should be exercised in 
the quality of the grid, the accuracy of inlet conditions, and the selection of the turbulence 
model. 

The analysis of the aerothermal field due to a compound angle leading edge film 
cooling indicates the presence of the complex vortex structure much different from that 
observed for a flat plate. Interaction between the upper coolant jet and the mainflow 
generates four major vortices; fXt a , fltb , & 5 » > £2 5b . These vortices originate respectively 
from the interaction of the mainflow with the coolant flow eminating from the right lateral 
side, core, left lateral side, and the upper side of the hole (Fig. 7-9 - Fig. 7-12). The 
vortices fXjb and £2 5b decay rapidly and the vortex flj b changes its sign beyond s/d=l 1. 
Vortices Qs a and fXjj, are major contributors to aerodynamic losses and for a decrease in 
the adiabatic effectiveness through the entrainment of the hot fluid. 
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Fig. 7-4 Grid structure and boundary conditions 
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Fig. 7-9 Vortex structure due to the upper hole jet-mainstream interaction, a) 
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Fig. 7-10 Vortex structure due to the upper hole jet-mainstream interaction, b) 
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Fig. 7-1 1 Vortex structure due to the upper hole jet-mainstream interaction, c) 
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Fig. 7-12 Vortex structure due to the upper hole jet-mainstream interaction, d) 
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Fig. 7-17 Velocity distribution , Tu=0.5%, contours represent distribution of 
normalized temperature, h 
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Fig. 7-19 Laterally averaged adiabatic effectiveness, influence of the Mach number 



Fig. 7-20 Surface pressure distribution at z/d=5 



Fig. 7-21 Laterally averaged adiabatic effectiveness, influence of high turbulence intensity 
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a) Tu=10%, 



z/d 


Fig. 7-24 Surface adiabatic effectiveness 



separation bubble 


Fig. 7-25 Streamlines from the bottom hole, Tu=0.5%; 



Fig. 7-26 Streamlines from the bottom hole, Tu=10% 




Chapter 8 


THREE DIMENSIONAL FLOW IN TURBINE ROTOR 

Numerical investigation of the flow field in the Penn State rotor based on the 
utilization of the Navier-Stokes solver has been carried out to gain a better understanding 
of the secondary and the tip leakage flow development. Pressure gradient across the blade 
tip clearance results in the development of the jet-like tip clearance flow. Its interaction 
with the meanflow leads to the development of the tip vortex and modification of the 
casing secondary flow and vortex. Similar to the jet-flow interaction analyzed in the 
previous chapter, the vortex interaction and mixing result in additional losses. One of the 
objectives of this simulation is to identify features affecting the development of secondary 
and tip leakage vortices, as well as sources of secondary flow losses. 

8.1 Computational details 

The current investigation is an extension of the numerical simulation presented in 
Luo and Lakshminarayana (1997). The emphasis of this research is to improve the 
resolution of the flow in the tip region, including a detailed analysis of the tip vortex 
development and assessment of the utilization of an ARS turbulence model. 

The embedded h-grid is utilized for the simulation of the flow in the tip region. 
There is a significant variation of velocity vector across the tip clearance due to the 
relative motion of the blade and casing. The maximum change in velocity amplitude 
reaches 90 m/s within 1mm distance. Thus, very dense grid is required for an accurate 
resolution of the tip clearance flow. A comprehensive grid dependency analysis is not 
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feasible due to the enormous utilization of the CPU time. Only a partial grid dependency 
analysis has been carried out. The numerical simulation has been performed with 9, 12, 
and 18 grid points in the tip clearance height (0.9 mm, T c =(rcas-rtip)/(r C as-rhub)=0.75). The 
difference between the solution based on 12 and 18 tip clearance grid points is 
considerably smaller than that based on 9 and 12 grid points within the gap. All results 
reported are based on the grid with 18 grid points in the gap. Similar analysis has been 
performed for the vortex region. Twenty- six grid points in spanwise-direction are utilized 
across the tip vortex zone. The total grid size is 104 (axial) X 60 (blade-to-blade) X 78 
(radial). The outlet boundary is set at x/C x = 2 downstream of the leading edge to ensure 
that there is no influence of the outlet boundary condition on the vortex structure 
development. Pitchwise average flow field based on the data by Zaccaria and 
Lakshminarayana (1995) is used to establish the inlet boundary conditions. 

8.2 Comparison with the experimental data 

Predicting strong secondary flow that exists in a turbine is known to be one of the 
most difficult tasks in CFD analysis. Results of the numerical modeling of the 
ERCOFTAC turbine test cascade reveal a large variation in the position and the amplitude 
of the secondary vortex predicted by various codes (Gregory-Smith, 1997). Current 
prediction of the flow in the Penn State rotor have been compared with LDV and pressure 
measurements by Ristic et al. (1998) and Xiao (1999) to establish confidence in the 
simulated results. 

A previous numerical simulation of the flow in the Penn State turbine Rotor (Luo 
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and Lakshminarayana, 1997) showed very good correlation between the predicted surface 
pressure distribution and design values from Hb=0.13 to H b =0.90. Comparison presented 
by Xiao, 1999, also indicates very good correlation between the experiment and predicted 
values. The numerical simulations presented in this thesis, which are based on a refined 
grid, variable tip gap, and a k-e /ARSM model, show that all these factors do not affect 
surface pressure at the blade location away from the casing region Hb<0.94. Thus that 
comparison is not presented here. However, the pressure distribution in the vicinity of the 
blade tip is strongly affected. Grid refinement has the most profound effect on the 
predicted pressure distribution near the tip clearance on the pressure side (Fig. 8-2). 

The comparison presented in Fig. 8-2 shows very good correlation between the 
predicted and measured blade pressure distribution on the pressure side. The presence of 
the tip clearance flow affects the pressure distribution only from 97% of span to the 
casing. The pressure field is essentially two-dimensional in nature (i.e., no significant 
variation in radial direction) from H b =50% to 97%. From 97% of the blade span, the 
velocity field undergoes sudden turning and acceleration as the flow enters the tip 
clearance. This leads to a rapid decrease in blade pressure near the tip from x/C x = 0.4 to 
x/C x = 0.9 in this region. This zone corresponds to the maximum tip-leakage mass flow at 
the pressure side of the gap. The low-pressure zone is relatively thin and extends from 
Zg=(r-r tip )/(r cas -rup)=-0.2 to Zg = 0.5. Grid density should be adequate not only in tip zone, 
but also below the blade tip to provide an accurate prediction of the flow in the flow 
acceleration zone. From Zg=0.5 to Zg=l (casing) the pressure level is relaxing to the 
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levels on the blade pressure surface. This is caused by the relative motion of the blade and 
casing, which generate the blockage near the casing and prevents or reduces the tip flow 
acceleration. 

Predicted pressure distribution on the suction surface also correlates well with the 
experimental data (Fig. 8-3). The characteristic feature of the suction side pressure field is 
the development of the low-pressure zone at x/C x ~0.5 near the blade tip. This pressure 
decrease is caused by the development of the tip vortex. The location of the minimum 
pressure can be used to identify the trajectory of the tip vortex (e.g., Ho and 
Lakshminarayana, 1996). A comparison between the predicted and measured location of 
the tip vortex path reveals earlier initiation of the tip vortex development in the case of the 
numerical simulation (x/C x ~0.58 VS x/C x ~0.64). At the blade trailing edge, the predicted 
distance between the vortex location and the endwall is higher than the measured value by 
about 1.5-2% of the blade span. One of the potential factors contributing to this 
discrepancy is the variation of the tip gap in the experiment. Design tip clearance is 
T c =l.l%, while the actual clearance varies from t c = 0.61% to t c = 0.9%. Base calculations 
have been carried out with -t c =0.75%. Additional simulations discussed below have been 
performed to assess the influence of the small variation of the tip gap height on the 
predicted flow field. 

Distribution of the axial velocity and axial vorticity at 10% of the chord 
downstream of the trailing edge is shown in Fig. 8-4 and Fig. 8-5. A comparison of the 
velocity field with the LDV data (Ristic et al., 1998) shows that the numerical solver 



232 


correctly predicts reduced axial velocity zone near the suction side caused by the 
secondary and leakage flows. The zone of the reduced axial velocity extends to 50% of 
the pitch near casing and to 80% of the pitch near the hub. Both, the experimental data 
and the prediction reveal the zone of the reverse flow in the center of the tip vortex (Fig. 
8-4). Distribution of the axial vorticity field shown in Fig. 8-5 may be used to identify 
major vortices; 1) casing wall passage vortex, 2) hub wall passage vortex, 3) tip leakage 
vortex, 4) scraping vortex (i.e. vortex caused by the relative motion of blade tip and 
casing) and 5) wake axial vortices. High flow turning results in the development of the 
passage vortices of significant strength. The hub wall passage vortex core is located at 
H b = 0.3 and spreads one third of pitch in tangential direction. LDV data has two zones of 
the positive vorticity near the casing; la and lb. A similar distribution can be observed in 
the predicted flow field (Fig. 8-5). However, measured axial vorticity in Zone la is 
significantly higher in comparison with the axial intensity in Zone lb. Predicted axial 


vorticity distribution has an opposite trend. Normalized helicity, 


mw 



can be utilized to 


analyze vortex development (e.g., Kunz and Lakshminarayana, 1992). Normalized helicity 
tends to unity at the vortex center disregarding the vortex intensity. The distribution of the 
predicted normalized helicity indicates that at x/C x ~0.8 the core of the casing wall passage 
vortex (defined as region with normalized helicity equal to unity) is separating into two 
parts. One is located closer to the casing and can be tracked to Zone la. The other part is 
transported downward and can be tracked to Zone lb. Based on this analysis, both la and 
lb are considered as elements of the casing wall passage vortex. A comparison between 
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the predicted flow field in the rotor with tip clearance and without tip clearance indicates 
that this separation of the casing passage vortex into two zones exists even in the absence 
of the leakage flow. Both, the prediction and the experimental data, clearly show that the 
tip vortex is confined to the suction side tip comer. The maximum leakage vorticity occurs 
at H b = 92% (prediction) and 94% (experiment). This correlates with the zone of low 
pressure on the suction side of the blade, which is discussed above. The tip clearance 
vortex extends to 20% in the pitchwise direction. The numerical simulation predicts a 
narrower tip clearance vortex. The wake development downstream of the trailing edge is 
three-dimensional in nature with a negative axial vorticity from 60 to 80% of the span and 
a positive axial vorticity in the lower 45% of the span. Interaction between the secondary 
flow and the wake, augmented by the rotation effects, is the primary mechanism 
responsible for the axial vorticity generation in the wake. Predicted wake axial vortices are 
more narrow than those observed in the experiment. This can be attributed to the 
interaction of the upstream nozzle wake with the rotor wake. This interaction, discussed in 
Chapter 5, makes the suction side of the rotor wake significantly thicker. The current 
simulation assumes a steady flow, thus the rotor wake-nozzle wake interaction effect is 
not captured. 

A comparison between the predicted flow field with the preliminary results of the 
LDV measurements (Xiao, 1999) is shown in Fig. 8-6-Fig. 8-9. The axial vorticity field is 
mostly two-dimensional distribution at x/C x = 50%, with the exception of the narrow zone 
near the casing, which corresponds to the development of the casing wall passage vortex. 
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Presence of a small zone of negative vorticity close to the blade tip indicates the inception 
of the tip clearance vortex. A comparison of the surface pressure distribution, presented 
earlier, indicates that the numerical simulation predicts earlier inception of the tip leakage 
vortex in comparison with the experiment. Therefore, predicted axial vorticity field at 
x/C x =50% of the chord has a larger zone of high axial vorticity. At x/C x = 80% (Fig. 8-9), 
the tip leakage vortex has grown significantly. The predicted axial velocity field indicates 
the presence of a stagnation zone in the core of the tip leakage vortex. This effect is 
weaker in the experimental data. However, the experimental data downstream of the 
leading edge (Fig. 8-4) contains a significant zone of the negative axial velocity, which is 
similar to the predicted flow. This discrepancy can be explained through a consideration of 
the tip vortex development. Presence of the reverse flow in the center of the tip leakage 
vortex, in combination with the fact that until 80% of the chord most of tip leakage flow is 
not entrained by the tip vortex, may lead to fewer LDV seeding particles in the center of 
the tip vortex. The tip leakage flow is essentially unsteady, therefore an experimental error 
may occur due to the variation of the position of the vortex core. 

An overall comparison between the predicted flow field and the experimental data 
is good and enables a certain level of confidence in the predicted flow field required for 
the analysis of the secondary flow development presented below. The discrepancy noted 
may be attributed to two factors. First is the limitation of the flow model utilized, i.e., the 
steady state simulation based on the circumferential average of the rotor inlet flow field. 
Experimental measurements (Lakshminarayana et al., 1998) and two-dimensional unsteady 
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simulation (Chapter 5) reveal a profound effect of the unsteady interaction on the flow 
field in rotor. Therefore, a steady state assumption may not capture some of the time- 
averaged features of the flow. The second factor is associated with deficiencies of the 
turbulence model. Ristic et al. (1998) measured strong anisotropy of the turbulence field 
downstream of the trailing edge. A further analysis is required to analyze the potential 
influence of non-isotropic turbulence on the accuracy of the prediction. 

8.3 Vortex field development 

Endwall boundary layers upstream of the rotor undergo strong modification as 
they enter the blade passage. Low momentum fluid located closer to the wall is 
transported towards the suction surface due to the blade-to-blade pressure gradient and 
streamline curvature. After the flow reaches the suction side, it develops into a passage 
vortices. Streamline and streamwise vorticity distribution (Fig. 8-10) illustrates the 
development of the secondary vortices in the rotor passage. Hub wall secondary flow 
impinges on the suction side of the blade at about 50% of the chord. At this location, it 
merges with the weak suction side horse shoe vortex. Further downstream, the hub wall 
passage vortex is transported away from the hub wall. At the trailing edge crossplane, the 
center of the hub vortex is located at 35% of the span. Development of the casing wall 
passage vortex is affected by the blade casing relative motion and by the leakage flow. 
Relative motion of the tip endwall contributes to a more intense transport of the casing 
boundary layer to the suction surface, in comparison with cascade flow. Casing boundary 
layer starts its final transformation into the passage vortex at x/C x = 40% of the chord near 
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the suction surface (Fig. 8-10). Starting at x/Cx= 50% of the chord, the development of the 
tip leakage vortex affects the casing wall passage vortex intensity and pushes it away from 
the suction surface. 

Relative motion of the casing wall results in a flow blockage (in pitchwise 
direction), preventing a development of the leakage flow within the first 30% of the chord 
and this strongly reduces the leakage flow further downstream. At x/C x = 20%, the 
tangential transport of the casing boundary layer fluid results in a reverse leakage flow 
from the suction surface to the pressure surface (Fig. 8-11, streamlines originated at 
x/C x = 20% of the chord and Fig. 8-13). The fluid, located closer to the blade tip, is 
transported along the blade centerline and leaves the tip gap only at x/C x = 50% of the 
chord (Fig. 8-14). This observation is supported by the distribution of the accumulated 
massflow rate through the gap plotted in Fig. 8-12. There is a weak negative tip leakage 
massflow rate up to 30% of the chord. Vector field at Zg = 0.33 indicates that the leakage 
flow entering the gap at the pressure surface near the leading edge leaves clearance only at 
x/C x = 50%. Within 60% of the blade, a significant part of tip leakage flow originates at 
the suction surface rather than at the pressure side of the blade tip as a result of the casing 
wall crossflow boundary layer. Near the casing, this zone extends from the leading edge to 
x/C x = 60% (Fig. 8-13) of the tip suction side. It shrinks to about 20% of the chord, from 
x/C x =40% near the blade tip (Fig. 8-14). Further downstream, an increasing pressure 
gradient across the tip gap confines a zone of the reverse leakage flow to less than 5% of 
the tip clearance height near the casing wall. At x/C x = 55% of the chord, all streamlines 
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initiated inside the gap propagates to the suction side of the blade, an indication of the 
normal pattern of the tip leakage flow (Fig. 8-11). The tip clearance flow originating at 
this location mixes with the mainflow without rolling up into a vortex. 

A strong interaction between the leakage flow and the casing wall boundary layer 
immediately turns the leakage flow downwards and streamwise as it leaves the clearance 
at x/C x = 60% (Fig. 8-14 and Fig. 8- 16a). An analysis of the secondary velocity field 
indicates the presence of a vortex inception between the blade surface and the leakage 
flow. Further downstream, this vortex is transformed into a full scale tip leakage vortex 
(Figs. 8-16). However, streamlines paths clearly show the absence of the tip leakage fluid 
inside this vortex at x/C x = 60% (Fig. 8-11 and Fig. 8-16). 

The inception of the tip leakage vortex occurs around 50% of the chord, caused by 
the interaction between the tip leakage flow and the main flow. Increasing leakage 
massflow expands tip leakage jet penetration into the main flow (Fig. 8-16). However, as 
stated above, the leakage flow does not start to roll up into the vortex until 80% of the 
chord. An enlargement of the crossflow area between the leakage jet and the blade suction 
surface leads to a significant de-acceleration of the flow in this zone. Ultimately, the zone 
of weak reverse flow develops at the center of the vortex zone. Streamlines initiated at the 
location of the flow separation zone shown in Fig. 8-17 indicate that the core of the tip 
vortex comprises of the mainflow fluid rather than the tip leakage fluid until 85% of the 
chord. The stagnation zone in the core of the tip leakage vortex grows steadily from its 
initiation up to the trailing edge (Fig. 8-16). Downstream of the trailing edge, the zone of 
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the reverse velocity disappears rapidly as a result of intense entrainment of the tip leakage 
flow into the core of the tip leakage vortex. 

From 80% of the blade chord to the trailing edge, there is a continuous change in 
the amount of the leakage flow entrained by the tip leakage vortex and the tip leakage 
flow mixing with the mainflow as a plain jet. At x/C x = 80% of the chord, only 10% of the 
tip leakage massflow ends in the tip leakage vortex (Fig. 8-16 and Fig. 8-18). At this 
location, most of the tip leakage fluid undergoes a quarter rotation as the outer layer of 
the tip leakage vortex and is then pushed downwards into the mainflow. Near the trailing 
edge, this part of the leakage flow is merging with the rotor wake (Fig. 8-18). Within the 
last 7% of the blade chord, practically all the tip leakage flow rolls up into the tip leakage 
vortex. The only exception is a thin zone near the casing wall. In this region, the scraping 
vortex prevents the flow from merging with the tip leakage vortex. 

Accumulated massflow through the gap grows rapidly from 40% of the chord (Fig. 
8-12). Maximum massflow rate is achieved from x/C x = 70% up to the 90% of the chord. 
The percentage of the leakage flow, which is entrained by the tip vortex, increases nearly 
linearly from x/C x = 80% of the chord to the trailing edge. Thus, the combined leakage 
massflow is equally split between the fluid, which is entrained by the tip vortex and the 
fluid which interacts with the mainflow as a plain jet. 

8.4 Secondary and leakage flow losses 

The complex structure of the vortices in the rotor and their interaction has a major 
influence on flow losses. There are three major sources of three-dimensional losses. The 
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first source of losses is due to the presence of strong casing and hub passage secondary 
vortices. The second contributor is additional losses associated with the development of 
the tip leakage flow. The last source is the increased or, in some cases, decreased losses 
due to the interaction between the tip leakage vortex, secondary vortices, and the rotor 
wake. 

The various zones of loss generation due to the leakage flow can be classified as: 

1 . Loss generation inside the tip gap, including losses due to the sudden contraction of the 
flow, tip and casing boundary layers, and a potential development of the separation zone; 

2. Mixing losses inside the blade passage. These losses occur due to the dissipation of the 
tip leakage vortex and “plain” leakage jet-mainflow mixing loss; 

3. Loss production associated with tip vortex development downstream of the trailing 
edge. 

The numerical simulation is a valuable tool in the investigation of the sources of 
additional losses, as well as their distribution. However, the predicted losses based on 
CFD modeling are not very reliable in terms of their absolute values. No experimental data 
is available at this time; therefore, predicted losses can not be verified. Nevertheless, the 
author’s experience and the information presented in literature show that CFD analysis 
provides reliable information in predicting the trend in loss distribution. 

Axial distribution of the mass averaged loss coefficient, 
(£ = (^>i -P 0 )/(p,Wj 2 /2), is presented in Fig. 8-19. A comparison with the loss 
coefficient based on the two-dimensional simulation, presented in Chapter 5, shows that 
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the secondary flow and tip leakage losses are responsible for about 50% of total rotor 
losses. This conclusion correlates well with experimental observations (e.g., Booth, 1985). 

The presence of the secondary flow and the development of the tip leakage flow 
(Fig. 8-20) results in increased losses downstream of the trailing edge, while the presence 
of the axial vortices in the wake contributes to the increased level of losses in this zone. 
This observation is similar to those made by Ho and Lakshminarayana (1996), for the 
turbine cascade. Inside the blade passage, increased rate of loss generation, in comparison 
with two-dimensional flow, is observed from 50% of the chord as a result of the final 
entrainment of the casing and hub wall secondary boundary layers into corresponding 
passage vortices. Intensified loss production from x/C x = 90% correlates well with the 
changing pattern of tip leakage - mainflow mixing. At this location, most of the leakage 
flow is entrained by the tip leakage vortex, resulting in additional losses. 

Losses inside the tip clearance gap 1 are responsible for about 6-7% of the total 
additional losses (Fig. 8-19). This low value may be attributed to the absence of the flow 
separation inside the gap and on the blade pressure surface in the vicinity of the gap. 
However, the prediction of the separation zone is sensitive to the characteristics of the 
turbulence model used. 

The secondary flow vortex structure is directly related to the distribution of the 
loss coefficient presented in Fig. 8-20. Development of tip and hub secondary vortices, 


1 Gap losses are calculated as the difference in the mass average stagnation pressure between the suction 
and the pressure sides of the gap, normalized by the ratio of tip leakage massflow to the total massflow 
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which can be observed at x/C x = 60% results in a zone of the low pressure near the tip and 
the hub wall respectively. Farther downstream, mixing of the secondary vortex is the 
primary source of losses up to the trailing edge. In the case of Penn State rotor, the 
presence of the tip leakage flow does not affect the loss generation due to the casing and 
hub passage vortices. This conclusion is based on a comparison of the loss coefficient 
distributions for the cases with and without tip clearance (Fig. 8-20 and Fig. 8-21). At 
x/C x = 80%, the loss coefficient distribution has a zone of the decreased total pressure 
located between the tip vortex center and the boundary layer (Fig. 8-20). Tip leakage flow 
streamlines initiated at x/C x at 70% of the chord inside the gap indicate that these losses 
are due to the mixing of the tip leakage jet (i.e., part of tip leakage flow mixing with main 
flow without entrainment into the tip leakage vortex). Beyond this zone there is no 
indication of the additional losses caused by the mixing of the leakage flow outside the tip 
leakage vortex. 

The development of the tip leakage vortex results in an extended zone of low 
pressure near the blade tip. The presence of the reverse flow in the vortex core minimizes 
the contribution of the zone of low total pressure inside the tip leakage vortex into overall 
losses. Downstream of the leading edge, massflow associated with the zone of the tip 
leakage vortex is higher. Therefore, the contribution of this zone in loss generation is more 
profound downstream of the trailing edge. 

Due to the complex interaction between secondary vortices, leakage flow, and 
wake, it is very difficult to calculate contribution of various sources to total loss. It is 
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estimated that the tip and the hub secondary vortices contribute 50% of the additional 
losses (i.e., losses in addition to the profile and plain wake mixing losses), while losses due 
to the tip leakage flow account for about 35% of the additional losses. The rest of the 
additional losses is due to the interaction of secondary vortex, the leakage flow and the 
wake (i.e., loss due to the presence of wake axial vortices, etc.). 

8.5 Influence of the tip clearance height on the tip leakage vortex 
development. 

A numerical simulation of the turbine rotor with different tip clearance heights has 
been carried out. The main objective of this investigation is to study the influence of 
clearance height on the amplitude and the structure of the secondary and leakage flows. In 
addition to gaining a better understanding, this simulation is useful in establishing a better 
interpretation of the comparison between the numerical and experimental results. This 
study will also provide guidance to the designer in optimizing the tip leakage effect. In the 
experimental rotor, the blade-to-blade variation of tip clearance is about 30% of its 
average value. Three cases have been investigated: with x c =l.l%, x c =0.61%, and no 
clearance (in addition to the base case with x c =0.75%). Even though zero clearance is 
physically not realizable, the relative motion of the blade and casing is preserved. A 
comparison of the leakage massflow at different tip clearances (Fig. 8-12) shows a 
nonlinear relation between the clearance height and the leakage massflow. The blade tip 
boundary layer and the casing boundary practically merge for the case of x c =0.61%. 
Therefore, leakage massflow at x c =0.61% is equal to less than one-third of those observed 
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at the design value of tip clearance (t c = 1.1%, Fig. 8-12). 

There is no significant difference (beyond the moderate change in amplitude of 
vortices and position of the tip leakage vortex core) between the axial velocity and axial 
vorticity fields at x/C x = 110% of the chord (Fig. 8-22 b,d and Fig. 8-4, Fig. 8-5) at 
various tip clearances. The position of the tip vortex center is located farther from the 
casing for the flow with x c =l.l%. This is due to the earlier inception of the tip leakage 
vortex caused by increased leakage flow. Based on the results of the numerical modeling, 
it is possible to estimate that 0.1% increase in clearance results in 1% change in the 
spanwise position of. the tip leakage vortex. Despite the significant change in the size and 
amplitude of the tip vortex, as well as in the extent of the tip leakage vortex core reverse 
flow, the overall pattern of the secondary and leakage flows in the case of x c =0.61% is 
similar to the base case (Fig. 8-22). However, there is a significant redistribution between 
the leakage flow entrained by the tip leakage vortex. For the case with t c =0.61%, about 
75% of the leakage flow is entrained into the vortex, while the base case has only 50% of 
the leakage flow rolled up into it. Mixing of the tip leakage vortex results in higher losses 
in comparison with “plain” jet mixing. 
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Fig. 8-1 Penn State rotor, computational grid 
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Fig. 8-3 Pressure distribution , suction surface (tip clearance zone is also shown in the 
prediction) 
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Fig. 8-4 Normalized axial velocity field at x/C x =1.10 
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Fig. 8-5 Normalized axial vorticity at x/C x =1.10 
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Fig. 8-18 Normalized streamline vorticity (oV2co), crossflow planes are located at 
x/C x = 60%, x/C x =80%, and x/C x = 1 10%, streamlines are initiated inside tip gap at 70% and 
90% of the chord 
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Fig. 8-21 Loss coefficient, no tip clearance 
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Fig. 8-22 a)-d) Flowfield at x/C x =l 10%, influence of tip clearance height a) Normalized 
vorticity, 0.61% clearance a) Normalized vorticity, 1.1% clearance, c) Normalized axial 
velocity, 0.61% clearance, d) Normalized axial velocity, 1.1% clearance 
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Fig 8-22 e)-f) Flowfield at x/C x =l 10%, influence of tip clearance height e) Normalized 
vorticity, no tip clearance, f) Normalized vorticity, no tip clearance 




Chapter 9 


SUMMARY 


9.1 Summary 

A compressible unsteady Navier-Stokes solver has been developed to enable 
numerical simulation of complex turbomachinery flows. To improve the computational 
efficiency of the code, a pseudo-time acceleration technique has been employed in the 
conjunction with the explicit Runge-Kutta scheme. Utilization of the pseudo-time 
approach leads to the presence of additional source terms, thus affecting both the stability 
and the convergence characteristics of the original scheme. Von Neuman analysis has been 
carried out to assess different methods of the pseudo-time stepping implementation. Based 
on this analysis, a correction to the local time step required for the stable and efficient 
unsteady calculations has been established. Quality of the grid plays a significant role in 
the accuracy of the numerical prediction. The ability to use different grid topologies, 
which are more suitable for the particular flow condition, is beneficial for improved 
prediction capabilities. It is also helpful in reducing user’s effort to generate the grid. The 
multiblock version of the code has been developed to make the code more flexible and 
more suitable for the numerical simulation of complex turbomachinery configurations such 
as multistage rotor-stator interaction, multidomain structure associated with film cooling. 


tip clearance flow etc. 



268 


Results of the numerical simulation based on the hybrid ARSM/ k-e model depend 
on the k-e component of the model, especially in the near wall region. This is especially 
noticeable in the transitional flows. Low Re number k-e model based on distance to the 
wall and wall shear stress was replaced with the k-e model based on local turbulence 
characteristic. This modification enables a better prediction of transitional and separated 
flows. Transition models have been incorporated in the code. This development is aimed 
at assessment of the numerical simulation of transition based on transition models against 
the prediction based on original turbulence models. 

The numerical solver has been extensively validated against benchmark flows. A 
number of criteria, required for the accurate numerical simulation, have been established 
based on the results of the validation. A comparison of the predicted turbulence field with 
the experimental data indicates that the present solver strongly overpredicts the level of 
the turbulence intensity in regions with the dominant normal stress. The modification of k- 
6 model has been performed to reduce the predicted level of the turbulent kinetic energy 
near the leading edge and away from the blade surface. The Navier-Stokes procedure has 
been used to investigate the unsteady flow in turbine and compressor cascades. The nature 
of the upstream wake propagation through the stage, and its interaction with the blade 
passage flow have been analyzed. The contribution of various components (potential 
interaction, viscous dissipation, inviscid stretching) into the overall wake decay has been 
discussed. The assessment of low turbulence models for their ability to predict the 
development of the unsteady transitional boundary layer has been carried out. 
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Numerical modeling of the transition over a laminar separation bubble has been 
chosen to assess the ability of the numerical solver based on different turbulence models to 
accurately predict the complex flow phenomena associated with the flow in a LP turbine. 

Finally, two cases of the flow with the complex vortex structure have been 
simulated. Numerical simulations of the leading edge film cooling and the tip vortex flow 
are aimed at a better understanding of the complex vortex flow caused by jet-mainstream 
interaction. Results of the numerical investigation are used to derive a better 
understanding of flow physics. Vortex structure is analyzed to identify sources of the 
aerodynamic losses and degradation in the heat transfer efficiency 

9.2 Conclusions 

Incorporation of the pseudo time stepping enables a significant improvement in the 
code performance. Required CPU time is from 5 to 25 times less than the CPU time 
required for the basic code, depending on the case. Utilization of the pseudo-time 
stepping is especially beneficial for the simulation of wake-blade row interaction. Navier- 
Stokes simulation requires only two-to-three times more CPU time in comparison with the 
Euler simulation, if the surface pressure distribution is of primary concern. Thus, the 
Navier-Stokes solver can be used as a replacement for the Euler code in a coupled 
Euler/boundary layer procedure. This will combine the efficiency and the accuracy of the 
unsteady boundary layer code with the ability of the Navier-Stokes code to predict 
accurate pressure distribution, upstream wake decay and the off-design and separated 
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flow. The required CPU time in this case is about three times as much as the CPU time 
required for the original Euler/boundary layer solver. 

The main result from the validation process is the establishment of the acceptable 
level of the artificial dissipation coefficient for accurate and stable solution. To provide an 
accurate flow simulation kt should be in a range from 0.005 to 0.015. This value of k 4 (the 
forth-order dissipation coefficient) is also adequate to provide the numerical prediction of 
the transition flow which is independent of the artificial dissipation. To obtain an accurate 
simulation of the upstream wake propagation through the turbomachinery stage, the grid 
should have at least twenty grid points per wave. However, this requirement can be 
moderately relaxed for high harmonics of the narrow wake (i.e., near wake). Rapid 
physical decay results in faster decay higher harmonics of the wake upstream of the 
leading edge. Thus the effect of artificial dissipation is not significant. 

The modification of the k-E model to eliminate overprediction of the turbulent 
kinetic energy near the leading edge and away from the blade surface is essential for 
accurate numerical simulation of transitional flows. 

9.2.1 Unsteady transitional flow in compressor cascade 

The ability of the Navier-Stokes code to predict the unsteady transitional flow on a 
turbomachinery blade is assesed. The unsteady pressure and velocity fields are in good 
agreement with the experimental data and the prediction from the Euler/boundary layer 
approach. The numerical solver is able to capture major zones (wake induced transitional 
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strip, wake induced turbulent strip etc.) associated with the wake induced transition in a 
compressor cascade. 

Another significant step is the assessment of k-e turbulence models, including 
leading edge modifications. Best results are obtained using FLB model. The LB model 
predicts earlier inception of the transition and shorter transition length. Modification of the 
k-e model is found to be essential for the accurate prediction of the unsteady transitional 
flow in the compressor cascade. The CH model fails to predict the unsteady transitional 
flow. The predicted boundary layer is turbulent from the leading edge, even with the 
modification of the k-e model near the stagnation point. 

The predicted momentum thickness reveals excellent agreement with the 
experimental data and the Euler/boundary layer prediction. Similar to the boundary layer 
solution, Navier-Stokes solver predicts higher level of time-average momentum thickness 
in comparison with the steady state solution. This is an indication of an increased loss due 
to the unsteady interaction 

Interaction between the upstream wake and the stator wake results in shedding of 
unsteady vortices from the trailing edge and increased dissipation in the stator wake and, 
as a consequence, increased rate of decay of the stator wake. 

9.2.2 Rotor-stator interaction in turbine stage 

Comparison of the predicted unsteady flow filed with the LDV data and dynamic 
pressure measurements indicates that the predicted velocity and pressure fields in a turbine 
stage are in good agreement with the experimental data at the design and the off-design 
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conditions. The measured unresolved unsteadiness and the predicted turbulence intensity 
show that the peak intensities are predicted reasonably well, but the wake width based on 
the unsteadiness shows that the computation has a larger diffusion (into the freestream) 
compared with the experimental data 

The pressure gust has the most influence near the leading edge. Its influence in the 
development of the unsteady flow-field is limited to 15% of the blade chord from the 
leading edge. Beyond this, the wake-blade interaction through the development of the 
counter-rotating vortices in the blade passage, is the source of unsteadiness. The 
maximum variation in the unsteady pressure was observed at x/Cx = 0.28 on the suction 
surface. 

Up to 15% of the chord, the viscous dissipation is responsible for 45% to 75% of 
the wake decay. Further downstream, the wake undergoes both the inviscid decay and the 
amplification inside the passage. The contribution of viscous dissipation is equal to 75% of 
the total decay and 58% of the decay in the passage. Most of additional losses due to 
unsteady flow occur upstream of the leading edge (~55%). Inside the blade passage, the 
nozzle wake mixing losses are small in comparison with blade profile losses. Increase in 
profile losses is attributed to the modification of the suction surface boundary layer due to 
the nozzle wake-blade interaction. 

The numerical solver was able to predict most of the features associated with the 
wake induced unsteady transition (wake induced transitional strip, turbulent strip, etc.). 
Even though the k-e model lacks physics of the bypass transition, the predicted 
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development of various zones follows the trend observed in the experiment. The major 
exception is the calmed region. The predicted region downstream of the wake induced 
transitional strip possesses some characteristics of the calmed region. However, this region 
can not be identified as a calmed region because it lacks other essential features of the 
calmed region, such as propogation of the zone trailing edge at O. 3 W 5 . 

The nozzle wake interaction with the rotor wake leads to an increased 
unsteadiness observed in both the experiment and the prediction. The suction side segment 
of the nozzle wake merges with the rotor wake, causing fluctuations in the rotor wake. 
The phase lag between the suction side segment and the pressure side segment of the 
nozzle wake gives rise to fluctuation in the rotor wake at double the nozzle wake 
frequency, but with a significantly smaller amplitude. 

There is a significant change in the overall flow pattern at the off-design condition. 
Due to an increased reduced frequency, there are more nozzle wake segments in the rotor 
passage in comparison with the design case. Another significant factor is the development 
of a strong separation bubble on the pressure surface within 15% of the chord near the 
leading edge. The interaction between the nozzle wake and the separation bubble results in 
amplified unsteady pressure on the pressure surface. Predicted unsteady velocity indicates 
that the nozzle wake-separation bubble interaction generates a vortex pattern in the 
separation zone, thus increasing unsteadiness in this zone. 
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9.2.3 Transition to turbulence over laminar separation bubble in LP turbine 

A numerical simulation of the flow in a low-pressure turbine was carried out to 
assess the ability of different turbulence models to predict transitional flow at different Re 
and turbulence levels. Best results have been obtained with the FLB model (without the 
transition model). Implementation of the hybrid k-e/ARSM improves the prediction for 
Re = 50,000, Tu = 2.5%. While having minimum impact in all other cases, this 
modification contributed to the redistribution of the turbulent kinetic energy between 
various components in the transition region. 

Utilization of the transition model does not result in an improved flow simulation. 
Analysis of the turbulence characteristics in the transition zone shows that the lack of 
improvement is due to interference between the transition model and the low Re 
turbulence model. In the current prediction, the transition inception from the “pure” k-e 
model is located only about 2% of the chord upstream of the measured location. An 
enforcement of the transition through the intermittency function leads to a double damping 
of turbulence in the transition zone. 

A number of factors have been found to be essential for an accurate prediction of 
the transition. The first factor is the need to limit the turbulence production near the 
leading edge to ensure an accurate development of the laminar boundary layer. 
Implementation of the fourth-order artificial dissipation in ' k-e equations, without 
modification for the leading edge flow, may lead to the development of the pseudo- 
turbulent boundary layer. The solver with a mixed second/fourth-order artificial dissipation 
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is less sensitive to this problem due to its ability to avoid an immediate transition inception 
near the leading edge at high levels of turbulence intensity. A second factor is the need to 
modify the freestream turbulence equation. This problem, as well as the first one, is due to 
the poor performance of a standard k-e model in the case of strong normal stresses. 
Without the adjustment of the freestream turbulence, the turbulence intensity may be over- 
predicted by 2-3%. The elevated level of the turbulent kinetic energy may affect the 
transition inception prediction. 

The establishment of the reliability range for the numerical solver is needed for its 
wider acceptance for the design problems. Even though grid independency has been 
verified through the numerical modeling, a further analysis shows that the predicted 
transition location is affected by the level of artificial dissipation. For small values of k 2 k e , 
the variation of the artificial dissipation acts as a disturbance (i.e., affecting the transition 
inception without diffusing k in the transition zone). This makes the assessment of the 
reliability of the prediction more complex. In addition to grid and turbulence model 
characteristics, the numerical scheme (i.e., form of differential approximation) and 
numerical details (e.g., the way the realizability of k is ensured) do contribute to the 
variation in the transition prediction. For the current solver, the potential error associated 
with this phenomenon can be estimated to be 2% of the chord for the location of inception 
of transition. For the attached flow, the level of accuracy achieved can be considered as 
acceptable even without utilization of the transition model. To ensure a reliable prediction 
of the separated flow cases, this level should be improved because even a moderate 
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variation in the predicted transition inception results in a significant variation of the 
separation bubble. 

9.2.4 Leading edge film cooling 

The numerical investigation indicates that the CFD analysis can be successfully 
employed for the prediction, simulation, and parametric study of the complex flows 
associated with the leading edge film cooling. However, great care should be exercised in 
the quality of the grid, the accuracy of inlet conditions, and the selection of the turbulence 
model. 

The analysis of the aerothermal field due to a compound angle leading edge film 
cooling indicates the presence of complex vortex structure much different from that 
observed for a flat plate. Interaction between the upper coolant jet and the mainflow 
generates four major vortices, pair of core vortices and pair of outer vortices. These 
vortices originate respectively from the interaction of the mainflow with the coolant flow 
emanating from the right lateral side, core, left lateral side, and the upper side of the hole. 
The outer vortices decay rapidly. Core vortices are major contributors to the aerodynamic 
losses and for a decrease in the adiabatic effectiveness through the entrainment of the hot 
fluid. 

Following conclusions can be made: 

1) The compound angle injection and the resulting large and differing flow 

angles of the jet and the mainstream result in large velocity gradient and a “jet-wake” 
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structure, leading to intense mixing. Large crossflow development results in the lifting 
of the jets and movement of hot spots away from the wall. 

2) Numerical simulation at a higher inlet Mach number shows only a small 
decrease in the cooling effectiveness due to the compressibility and modified pressure 
gradient effects. 

3) High turbulence intensity leads to a decrease in the cooling effectiveness 
near the stagnation surface. The development of the flow above the second row of 
holes at the high turbulence level is controlled by two opposing trends; increased 
turbulence dissipation and modification of the vortex structure leading to a decrease in 
the hot fluid entrainment. The freestream turbulent length scale has a significant effect 
on the balance between these two phenomenon. 

4) The predicted cooling effectiveness and the flow-field are very sensitive to 
the inlet turbulent length scale. Increased length scale results in decreased effectiveness 
and faster decay of the vortices. Specification of the correct turbulence dissipation 
ratio or the length scale is one of the crucial elements for an accurate prediction of the 
leading edge film cooling effects. 

A more accurate prediction of the film cooling effects may be achieved through an 
improved turbulence model. The incorporation of the higher order anisotropic turbulence 
model and the modification of the constant turbulent Prandtl number assumption may 
improve the accuracy of the vortex simulation. Additional investigations should be carried 
out to assess the influence of the curvature effect on this type of flow. 
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9.2.5 Secondary flow in turbine including tip leakage flow 

Development of the tip leakage flow is characterized by significant velocity and 
pressure gradients that exist in the tip gap and its vicinity. According to the result of the 
grid refinement analysis, it is essential to have at least 15 grid points within the gap for an 
adequate numerical resolution of the flow. The sudden contraction of the flow generates a 
low-pressure zone below the blade tip on the pressure surface. Therefore, grid also should 
be highly clustered at least two tip gap heights below the blade tip. 

There is a good correlation between the predicted and the measured blade pressure 
distribution. The comparison between the predicted flow-field and the LDV data also 
reveals good correlation downstream of the trailing edge. However, the numerical 
prediction indicates an earlier development of the tip leakage vortex. Pitchwise width of 
the tip leakage vortex is smaller in comparison with the experiment. These discrepancies 
can be attributed to the limitation of the physical model, especially to the steady state 
approximation and isotropic nature of the turbulent model. 

Development of the leakage flow in the rotor significantly differs from that 
observed in a cascade. Relative motion of the blade and casing blocks the development of 
the tip leakage flow. Massflow through the tip gap reaches its maximum at 85% of the 
chord. Leakage flow leaving the tip clearance is only partially rolled up into the tip leakage 
vortex. At 50% of the chord, all the leakage flow mixes with the mainflow as a “plain” jet, 
while at 93% all leakage flow leaving the gap is entrained by the vortex. Tip vortex starts 
at around 50% of the chord as a radial separation zone on the suction surface. It 
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immediately moves away from the blade surface and grows steadily until the trailing edge. 
The core of the tip leakage vortex mostly consists of the mainflow fluid entrained at the 
vortex inception up to the 90% of the blade. In the case of the intense tip leakage vortex 
growth, a zone of reverse flow develops in the core of the tip leakage vortex, increasing 
the flow blockage due to the tip leakage flow. 

A comparison between the three-dimensional prediction, and two-dimensional 
prediction of pitchwise mass average loss coefficient shows that the secondary and leakage 
flow losses are responsible for about 50% of total losses. The prediction shows a relatively 
small contribution of the tip leakage flow (less then 30% of the additional losses). This can 
be attributed to the relatively small clearance (less then 1% of the blade span) in Penn 
State rotor. Most of the losses due to the tip leakage flow occur downstream of the 
trailing edge through the tip vortex mixing. 

A smaller tip clearance results in decreased leakage flow. It also leads to an 
increased percentage of the tip leakage fluid entrained by the tip leakage vortex. 
Therefore, the decrease in losses is less prominent, because of higher losses associated 
with the flow entrined by the tip leakage vortex. 

9.3 Recommendations for the future work. 

Reliable prediction of the calmed region can be beneficial for the improved blade 
design, especially potential improvement in stage weight characteristics. Even though the 
accuracy of the unsteady transition flow prediction is comparable with the steady 
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transition prediction, some essential features of the calmed region are not captured. This 
includes such characteristics as the propagation velocity of the calmed region near the 
trailing edge. Improved turbulence modeling is needed in order to get better resolution of 
the wake induced transition. 

The numerical prediction of the separated flow in LP turbine is very sensitive to 
the accuracy of the transition inception. Small error in the predicted location results in 
considerably modified height and extent of the separation bubble. Further analysis is 
required to ensure a more stable and reliable simulation. Numerical solver should be 
modified to eliminate the dependency of the predicted transition on the numerical 
realization. 

Complex vortex structure due to the film cooling is known to possess non- 
isotropic turbulence nature. The ratio of eddy viscosity coefficients in streamwise and 
spanwise directions are significantly different from unity. This ratio can be as high as 20. 
Additional research is required to explore the effect of non-isotropy in turbulence on the 
accuracy of the vortex flow simulations. 

Present compressible numerical solver was successfully employed for the flow 
simulation with relatively low flow Mach number. However, code utilization for 
calculations of three dimensional, low speed flow results in code computational efficiency 
degradation, especially for cases with significant variations in the. total energy distribution. 
An incorporation of the preconditioning may be beneficial for computationally more 
efficient modeling of low Mach flows. 
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Appendix A 


LINEARIZED SOLUTION OF THE WAVE PROPAGATION 


Let us consider sinus wave propagation: 
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Governing laminar incompressible equations can be written: 
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Artificial dissipation (assuming isotropic dissipation) is given by 
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where u and v are x and y components of the total velocity. 
Solution for the eq. [A-l] can be found in the form: 
u = u 0 ( 1 + A(x)cos(2^0J(y - V 0w t)) 


v = V 0W - u 0 J_ C0S (2^7(y - V 0w t) [ A-3] 

dx 2m 

Solution for A(x) can be obtained by the substitution of [A-3] into [A-l], 

In this case equation for the x-component can be rewritten as: 
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where (•)=(y-V 0w t) 

Linearization assumption, A(x) «1, has been used to derive this expression. 

Cell aspect ratio is assumed equal to unity, i.e.. Ax ~ Ay 

Equation [A-4] can be non-dimensionalised using wave length :x = x/l , where 1 - is 
length of the wave. 
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Where Re w is Reynolds number based on the wave length. 
If CFL = 1 then 
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where n is the number of grid points per wave length. 

Making additional assumptions 

(27COJ) 2 A(x) 

OX 

^^«(27rGJ) 4 A(x) 

ox 

This assumption is based on the dominance of the crossflow component of the dissipation 
over a streamline one. 

Then, solution for [A-5] can be written as: 

, 1 (27TGJ) 2 

A(x) = 4 3 exp[(-2*G7) 2 (-— + k 4 ±— j+-)x] 

Re w n 

For the freestream wave/wake propagation, discussed in chapter 3, u and v are the 
functions of (y-V 0w t+VouAiox). This problem may be reduced to the problem considered 
above using coordinate transformation : t, = y — V 0w x 
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